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ABSTRACT 

A  Navier-Stok.es  problem  solver,  developed  by  L.  N.  Sankar,  is  modified  to  provide 
dynamic,  interactive  graphical  presentations  of  predicted  flow  field  solutions  about  a 
NACA-0012  airfoil  section,  oscillating  in  pitch,  in  order  to  demonstrate  the  capabilities 
of  dynamic  graphics  applications  in  the  study  of  complex,  unsteady  flows.  Flow  field 
solutions  in  the  form  of  pressure  coefficient  and  stream  function  contour  plots  about  an 
airfoil  experiencing  dynamic  stall  are  plotted  utilizing  an  IRIS  3000-series  workstation 
and  Graphical  Animation  System  (GAS)  software,  developed  by  Sterling  Software  for 
NASA.  These  full  cycle  solutions,  in  conjunction  with  dynamic  surface  pressure  dis- 
tribution plots  and  integrated  lift,  pitching  moment  and  drag  coefficient  data,  are  com- 
pared to  existing  experimental  data  in  order  to  provide  an  indication  of  the  validity  of 
the  code's  far-field  solution.  Full  procedural  documentation  is  maintained  in  order  to 
provide  an  efficient  analysis  tool  for  use  in  future  oscillating  airfoil  studies  planned  by 
the  NASA-Ames  Fluid  Mechanics  Laboratory  and  the  Naval  Postgraduate  School  De- 
partment of  Aeronautics  and  Astronautics. 
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I.     INTRODUCTION 

It  has  been  recognized  for  some  time  [Ref.  1]  that  an  airfoil  oscillating  in  pitch  to 
angles  of  attack  greater  than  the  static  stall  angle  will  surpass  the  traditional  stall  burner 
and  generate  normal  forces  which  exceed  those  attainable  in  the  static  ci.se.  This  dy- 
namic stall  phenomenon  is  attributed  to  the  aft  movement  of  a  strong  shed  voitex  along 
the  upper  surface  of  the  airfoil,  carrying  with  it  an  induced  velocity  Held  which  radically 
and  dynamically  changes  chordwise  pressure  distributions.  In  general,  an  accurate  in- 
terpretation of  the  dynamic  stall  mechanism  will  significantly  impact  a  variety  oi  appli- 
cations, all  of  which  involve  dynamic  lifting  surface  motions  and  unsteady  flow 
separation.  Originally,  dynamic  stall  analysis  efforts  were  directed  toward  helicopter 
aerodynamics,  where  sharp  increases  in  oscillatory  torsional  loading  and  thus,  blade 
stress,  can  reduce  the  fatigue  life  of  rotor  mechanical  components  and  vortex-induced 
aerodynamic  loading  can  generate  adversely  phased  pitching  moments,  resulting  in  stall 
flutter  [Ref.  2|.  Much  current  interest  concerns  the  feasibility  of  exploiting  dynamic  stall 
forces  for  effective,  sustained  maneuvering  in  the  high  angle  of  attack,  "post-stall"  flight 
regime,  or  supermaneuverability.  for  next-generation  fighter  and  attack  aircraft  [Ref.  3 J. 
Historically,  attempts  to  analyze  such  complex,  unsteady  behavior  relied  heavily  on 
empirical  data,  obtained  from  often  laborious  and  time-consuming  tests.  With  the  ad- 
vent of  the  supercomputer,  however,  computation  of  actual  viscous  flow  fields  about 
moderately  complex  computational  models  can  now  be  numerically  achieved  in  a  matter 
of  minutes  through  solution  of  the  Reynolds-averaged  Navier-Stokes  equations  [R.ef 
4].  These  solutions,  in  and  of  themselves,  however,  are  not  sufficient  to  promote  insight 
into  the  mechanics  and  physics  involved  in  such  flows.  This  additional  requirement,  for 
effective  visual  portrayal  of  the  flow  field,  is  satisfied  by  application  of  high  performance 
interactive  computer  graphics  workstations  and  associated  software. 

The  long  range  goal  of  the  Computational  Fluid  Dynamics  field  is  the  development 
of  a  thoroughly  verified  computer  code  for  unsteady  aerodynamics  which,  among  a 
myriad  of  other  applications,  will  provide  future  aircraft  designers  with  the  opportunity 
to  derive  full  advantage  from  application  of  the  dynamic  stall  phenomenon.  To  this  end. 
the  Fluid  Mechanics  Laboratory  (F.ML)  of  the  NASA-Ames  Research  Center  (ARC). 
in  conjunction  with  the  Naval  Postgraduate  School  Department  of  Aeronautics  and 


Astronautics  has  planned  an  assortment  of  parallel  and  complementary  oscillating  airfoil 
windtunnel  experiments  and  computer  simulations. 

The  current  study  utilizes  existing  dynamic  graphics  packages  to  generate  real-time 
projections  of  flow  fields  about  a  NACA-0012  airfoil  oscillating  in  pitch  and  experienc- 
ing deep  dynamic  stall,  as  provided  by  a  Navier-Stokes  solver  developed  by  L.  N.  Sankar 
[Refs.  5,6].  A  modified  version  of  the  Sankar  code  was  submitted  via  the  FML  front  end 
VAX,  to  the  NASA-ARC  Cray  X-MP/48  computer,  which  output  flow  field  solutions 
at  specified  intervals  throughout  the  oscillatory  cycle.  From  this  data,  graphics  files  were 
generated  which,  in  turn,  were  submitted  to  the  Graphics  Animation  System  (GAS) 
software  as  developed  by  Sterling  Software  under  contract  to  NASA- ARC,  and  run  on 
an  IRIS  3000-series  graphics  workstation.  (Final  output,  for  demonstrative  purposes, 
was  then  transferred  to  video  tape.)  Thus,  the  results  of  the  study  were  two-fold.  First, 
procedures  by  which  interactive  computer  graphics  could  be  efficiently  (in  terms  of  both 
computer-time  and  man-hours)  and  effectively  (in  terms  of  information  display)  incor- 
porated as  an  analysis  and  verification  tool  for  future  FML  studies,  were  developed, 
tested  and  refined.  Secondly,  the  resultant  graphics  were  utilized  as  a  tool  for  the  on- 
going verification  of  the  Sankar  code. 


II.     DESCRIPTION  OF  THE  SANKAR  CODE 

Simulation  of  complex  phenomena  occurring  in  real  fluid  flows  requires  accurate 
solutions  to  the  lull  Navier-Stokes  equations.  The  Sankar  Navicr-Stokes  solver,  devel- 
oped for  Blade-Vortex  Interaction  (BVI)  studies,  solves  the  two-dimensional,  unsteady, 
compressible  Euler  and  N'avier-Stokes  equations  in  strong  conservation  form,  utilizing 
an  alternating  direction,  implicit  method  as  the  time  marching  algorithm.  A  body-fitted 
C-arid,  with  clustering  in  the  normal  diiection  is  utilized  to  discretize  the  flow  field. 
Turbulent  shear  stresses  are  simulated  with  a  two-layer  algebraic  eddy-viscosity  model, 
with  modifications  as  described  later.  The  code  may  thus  be  utilized  for  solution  of 
steady  or  unsteady,  inviscid  or  viscous  and  laminar  or  turbulent  flows. 

A.     GOVERNING  EQUATIONS 

In  Cartesian  coordinates,  the  2-D,  unsteady,  compressible  Navier-Stokes  equations 
in  strone  conservative  form  mav  be  written  as 


5tq  +  SXE  +  SyF  =  Re   \oxR  +  5  S) 


(2.1) 


where  q,  E,  F.  R  and  S  are  four  element  vectors  with  entries  corresponding,  in  order,  to 
the  equations  of  continuity,  momentum  (x-  and  y-)  and  energy.  If  non-dimensionalized 
such  that  all  values  of  length,  velocity  (u  and  v),  density  (  p  )  and  total  energy  per  unit 
volume  (e)  are  normalized  with  respect  to  section  chord  length  (c),  free  stream  speed  of 
sound  (  aj),  frce  stream  density  (  px)  and  pxax,  respectively,  these  vectors  may  be  writ- 
ten as 


p 

P 

pv 

0 

0 

pu 

£  = 

pu2  +  p 

F  = 

puv 

2 

R  = 

* XX 

s  = 

^xy 

pv 

puv 

pv    +p 

*xy 

Tyy 

e 

_u{e+p)_ 

_v{e  +  p)_ 

,R*. 

s,_ 

(2.2) 


The  Reynolds  number  (Re),  pressure  (p),  speed  of  sound  (a)  and  stress  terms  (  tm)  are 
defined  as  follows. 


Re»  £  (2.3) 

;,  =  (y_  l)|>-0.5p(u2  +  v2)]  (2.4) 

a2  =  y(y  _  i)[ejp  -  o.5(u2  +  r2)]  (2.5) 

txx  =  {/.  4-  2,u)ux  +  ).\-y  (2.6) 

Txy  =  n(Uy  +  l\J 
Ty>,  =  (/  +  2/i)^  +  ).ux 
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Assuming  Stokes'  hypothesis  to  be  valid,  the  constant  i  is  defined  as  -2/3  \i  .  Since  only 
airflow  is  considered,  the  ratio  of  specific  heats  (  y  )  is  defined  as  1.4  and  the  Prandtl 
number  (Pr),  as  one. 

Neglecting  viscous  terms  results  in  the  Euler  equations,  as  the  right-hand  side  of 
equation  (2.1)  goes  to  zero.  L'se  of  the  strong  conservative  form  (i.e.,  the  continuity 
entry  in  E,  6{pu)jdx,  is  solved  in  its  present  form,  vice  the  more  mathematically  correct 
form  of  pdujdx  +  udpjdx)  allows  the  identical  conservation  of  physical  quantities  (mass, 
momentum  and  energy)  when  finite  difference  schemes  are  applied. 

B.     TRANSFORMED  GOVERNING  EQUATIONS 

The  flow  field  region  must  be  discretized  into  a  transformed,  finite  difference  mesh, 
or  computational  plane,  in  order  to  allow  numerical  solution  of  the  governing  equations. 
The  most  efficient  grids  are  rectangular  in  shape  with  regular  spacing  or  spacing  gradi- 
ents. They  must,  however,  correspond  to  a  flow  field  grid  which  provides  high  resolution 
(minimal  spacing)  in  regions  where  gradients  are  large,  particularly  in  the  boundary  layer 
and  about  the  leading  edge.  In  response  to  the  coordinate  transformation  involved  in 
the  grid  generation  process,  the  governing  equations,  too,  must  be  transformed.  By  de- 
fining £  and  t]  as  functions  of  the  cartesian  coordinates  (time  is  not  transformed),  this  is 
simply  a  2-D  mapping  procedure  accomplished  by  application  of  the  transformation 
Jacobian, 
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the  governing  equation  becomes 
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where 
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(2.10) 
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C.  GRID  GENERATION 

Grid  generation  for  discretization  of  the  flow  Held  region  may  be  accomplished 
conformal  mapping,  algebraic  methods  or  numerically  solving  a  set  of  partial  differential 
equations.  The  original  Sankar  code  utilizes  either  of  the  latter  two  methods,  while  the 
version  currently  vectorized  for  use  in  this  study  utilizes  only  the  the  algebraic  method, 
which  results  in  a  body-fitted,  sheared  parabolic  coordinate  system  or  C-grid.  Utilization 
ol  such  body-fitted  grids  allows  synchronous  rotation  of  the  entire  grid  and  an  foil  sec- 
tion for  dynamic  cases,  using  simple  trigonometric  relations.  This  coordinate  system 
satisfies  the  general  grid  requirements  for  smoothness  throughout  and  fine  spacing  in 
regions  where  high  gradients  exist,  such  as  the  boundary  layer  or  leading  edge. 

Normalized  geometric  airfoil  shape  data,  in  Cartesian  coordinates,  is  input  in  table 
format  to  the  code,  which  utilizes  an  interpolative  procedure  to  compute  additional 
points  and  smooth  the  surface.  The  trailing  edge  region  is  modelled  as  a  vortex  sheet 
shape,  or  "cut",  which  smoothly  leaves  the  airfoil,  tangent  to  the  mean  camber  line  at 
the  trailing  edge.  Algebraic  manipulation  of  the  section  surface  and  cut  allows  grid 
generation  in  the  transformed  £  -  ;/  plane,  figure  1  shows  the  resultant  grid  when 
mapped  back  to  Cartesian  coordinates. 

Uniform  spacing  in  the  transformed  plane's  c  -direction  results  in  fine  spacing  at  the 
leading  edge  in  the  real  plane,  but  coarse  spacing  in  the  nailing  edge  region,  due  to 
uniform  increments  across  the  axis  (as  opposed  to  the  standard  C-grid  which  results  in 
a  region  of  finer  spacing  at  the  trailing  edge).  >/  -spacing  in  both  planes  increases,  ap- 
proximately exponentially,  in  directions  normal  to  the  airfoil  surface,  with  the  magnitude 
of  initial  spacing  being  user  defined.  Though  easy  to  generate,  this  grid's  effectiveness 
in  the  analysis  of  certain  flow  cases  has  proven  somewhat  limited,  due  to  its  coarseness 
in  the  trailing  edge  region  [Ref.  6,  p. 44]. 

D.  APPLICATION  OF  BOUNDARY  CONDITIONS 

Consideration  of  an  airfoil  impulsively  started  from  rest  in  a  fluid  with  uniform 
properties  throughout,  provides  the  initial  conditions  for  solution  of  the  parabolic  Euler 
and  Navier-Stokes  equations  (Eq.  2.9)  in  the  computational  plane.  In  addition,  bound- 
ary conditions  and  artificial  dissipation  terms  are  required  in  order  to  achieve  accurate 
solutions. 

Three  boundaries  exist  which  are  of  interest,  the  surface,  the  far-field  boundary  (grid 
limit)  and  the  trailing  edge  cut.  Boundary  conditions  at  the  surface  are  driven  by  the 
no-slip  condition  which,  in  response  to  viscous  effects,  requires  the  fluid  velocity  at  the 
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Figure  1.      Algebraic  C-Grid  in  Cartesian  Coordinates 


boundary  to  equal  the  boundary  velocity.  Thus,  in  this  reference  frame,  u  and  v  are  set 
to  zero  on  the  solid  surface.  Since  the  surface  is  considered  adiabatic,  both  6e'6)j  and 
de/dC  are  set  to  zero  and.  likewise,  pressure  distributions  are  determined  from  Opi'dy  and 
dplS?  equal  to  zero  at  the  solid  boundary.  (This  is  equivalent  to  neglecting  stress  con- 
tributions in  the  momentum  equations,  which  is  acceptable  when  dealing  with  high 
Reynolds  numbers  flows.)  Boundary  conditions  at  the  cut  are  driven  by  the  necessity 
to  avoid  discontinuities  in  the  "continuous"  portion  of  the  How  field.  Since  the  grid  is 
extremely  dense  in  the  >i  -direction  in  this  region,  averages  of  the  values  of  the  two 
nearest  interior  points  are  assigned  to  points  along  the  cut,  allowing  a  smooth  transition. 
Since  the  grid  cannot  economically  be  generated  large  enough  to  avoid  disturbance  ve- 
locities at  the  far-field  boundary,  boundary  conditions  must  account  for  the  presence  oi' 
the  airfoil.  Linear  small  disturbance  theory  is  applied  to  determine  perturbation  veloci- 
ties at  points  along  the  boundary,  which  are  then  added  to  free  stream  conditions. 
Downstream  boundaries  are  treated  such  that  entropy  changes  can  be  convected  out  of 
the  computational  domain  [Ref.  6.  p.  29|,  allowing  shocks  and  boundary-layer  generated 
vorticity  to  pass  through  the  grid. 

It  has  been  found  that  spatial  derivatives  are  sensitive  to  the  decoupling  of  odd  and 
even  points  which  necessarily  occurs  in  central  difference  schemes.  This  results  in  the 
generation  of  high  frequency  errors  in  regions  of  large  pressure  gradients,  such  as  are 
present  in  shocks  or  about  stagnation  points.  Due  to  the  high  Reynolds  numbers  en- 
countered in  these  flows,  dissipation  provided  by  the  viscous  terms  are  not  enough  to 
eliminate  the  errors,  necessitating  the  addition  of  two  artificial  dissipation  terms,  em- 
bedded in  the  numerical  schemes.  An  explicit  artificial  viscosity  term  is  input  by  the 
user,  with  a  proportional  implicit  term  assigned  by  the  code. 

E.     TURBULENCE  MODELLING 

Since  during  the  derivation  of  the  Navier-Stokes  equations,  no  assumptions  are 
made  regarding  flow  type,  these  equations  are  instantaneously  valid  for  both  laminar 
and  turbulent  flows.  The  large  range  of  time  and  spatial  scales  encountered  in  turbulent 
flows,  however,  makes  solution  of  the  instantaneous  Navier-Stokes  equations  virtually 
impossible,  at  this  time.  As  a  result,  the  equations  are  time-  or  ensemble-averaged.  Use 
of  such  equations,  in  response  to  turbulence,  results  in  additional  turbulent  shear  stress 
terms  or  the  Reynolds  stresses  (named  for  Oswalde  Reynolds,  who  initiated  turbulence 
studies  in  the  1880's),  which  effectively  are  increases  in  shear  stresses  due  to  turbulent 
motion.    The  difficulties    associated  with  the  existence  of  additional  unknowns  without 


any  additional  equations  is  described  as  the  closure  problem  and  is  dealt  with  through 
the  use  of  turbulence  models.  These  are  based  on  empirical  knowledge  of  turbulent 
flows  and  thus,  provided  solutions  will  be  approximate  and  require  some  confirmation 
from  experiment  [Ref.  7|. 

The  most  common  turbulence  modelling  method  involves  mixing  lengths,  as  devel- 
oped by  Prandtl.  Assuming  that  turbulent  fluctuations  arc  essentially  the  result  ol~  ve- 
locity perturbations  between  adjacent  streamlines,  and  that  all  fluctuating  velocity 
components  at  a  given  point  are  of  the  same  magnitude,  it  can  be  shown  that  turbulent 
or  eddy  viscosity  (  nT  )  is  proportional  to  the  magnitude  of  the  local  vorticity  (  w  )  [Ref. 
8,  p.  387],  In  Prandrl's  mixing  length  model,  the  proportional  term  is  the  square  of  a 
characteristic  length,  related  to  fluid  turbulence  intensity.  Prandtl  suggested  this  char- 
acteristic length  (  /  )  be  treated  as 

l  =  ky 

where  y  is  the  normal  distance  from  the  fluid  boundary  and  k  is  empirically  obtained. 
Thus,  the  key  to  obtaining  accurate  solutions  when  utilizing  models  of  this  type  lies  in 
the  mixing  length  expression. 

The  Baldwin-Lomax  two-layer  eddy  viscosity  model,  based  on  Cebeci's  two-layer 
model,  is  presently  used  in  the  Sankar  code.  This  model  divides  the  boundary  layer  into 
two  regions,  an  inner  layer  and  an  outer  layer,  with  separate  methods  for  determining 
eddy  viscosity  used  in  each.  The  boundary  for  the  two  layers  is  defined  as  that  point 
where  eddy  viscosities  produced  by  the  two  methods  match.  1  he  inner  layer  utilizes  a 
mixing  length  model  where  eddy  viscosity  starts  from  zero  at  the  wall  and  is  defined  by 
the  following. 

(dinner  =  PM  w  I  t112^ 

The  mixing  length  term  is  defined  by 

1  =  kvD  (2.13) 

where y  is  the  normal  distance  from  the  wall,  k  is  the  Von  Karman  constant  and  D  is  the 
Van  Driest  damping  function,  given  by 

D  =  [l-exp(->-+M+)]  (2.14) 

where  A*  is  an  empirical  constant.   The  outer  layer  eddy  viscosity  model  is  defined  by 


(a*  i)outer  =  KCcppFwakeFK!ebtv)  (2.15) 

where  K  and  Ccp  are  constants.    The  K       inolTintermittency  function,  given  by 

FiaJy)  -  CI  +  5.5(Qw/>maxj6]-*  (2.16) 

ensures  that  eddy  viscosity  approaches  zero  as  the  edge  of  the  boundary  layer  is  ap- 
proached and  the  flow  assumes  external  characteristics.  CKUb  is  a  constant.  Fwakt  is  a 
function  defined  the  following  relation. 

Fwake  ~  min0;max'7max'  ^w/J'max^V^max)  (2- 17) 

Udi/  is  the  magnitude  of  the  velocity  profile's  velocity  range.  Fmix  is  the  maximum  value 
provided  by  the  following. 

F(y)=y\co\D  (2.18) 

ym3X  is  the  value  of  y  corresponding  to  Fmax.  Constants  in  the  Sankar  code  are  defined 
as  follows: 

A+  =  26.0 

k  =  0.4 

K=  0.0168 

Ccp=\.6 

Qw  -  0.3 
CwA  =  0.25 

In  order  to  account  for  turbulence  outside  the  boundary  layer,  such  as  that  which 
occurs  in  the  dynamic  stall  process,  a  modified  turbulence  model  is  available  which  ef- 
fectively increases  the  outer  model's  mixing  length.  In  this  case,  Fmtx  andj,'max  are  deter- 
mined by  the  following 

F(y)=y2\co\D  (2.19) 
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The  modified  turbulence  model  was  developed  in  response  to  early  predictions  of  flow- 
separation  at  high  angles  of  attack  and  has  been  found  to  cause  premature  reattachment 
during  the  downstroke  of  dynamic  cycles.  Its  use,  therefore,  is  only  recommended  until 
stall  onset  [Ref.  8. p.  33]. 

F.     CODE  STRUCTURE 

The  code  is  structured  such  that  the  following  user  defined  parameters  are  input 
from  logical  unit  05  (appended  to  the  end  of  the  code  for  Cray  processing):  grid  size,  step 
size  (dt),  artificial  viscosity  magnitude,  mean  angle  of  attack  (a0),  oscillation  magnitude 
(a,),  angle  for  suspension  of  the  modified  turbulence  model,  reduced  frequency  (k).  free 
stream  Mach  number  (M^),  Reynolds  number  (Re),  distance  of  the  first  r\  -  contour  from 
the  airfoil  wall,  starting  time,  pitch  and  restart  flags  and  airfoil  geometry  data.  Added 
to  this  list  for  the  present  study  are  the  number  of  time  steps  for  the  code  to  march  on 
the  present  run,  the  number  of  steps  between  each  plot  (output  interval)  and  the  total 
number  of  steps  and  plots  completed  on  all  previous  runs.  Variables  are  then  initialized, 
previous  stored  solution  datasets  are  retrieved  and  read  (to  allow  incorporation  of  all 
datasets,  including  those  from  the  present  run,  into  a  single,  combined  file)  and  flow  field 
starting  conditions  are  input.  A  series  of  subroutine  calls  then  generate  the  grid 
(AIRFOL,  SING,  TABINT,  WRAP),  cluster  grid  points  for  viscous  flows  (CLUSTR, 
STRTCH),  rotate  the  airfoil  (ROTGRID)  to  its  actual  angle  of  attack  and  compute  the 
initial  metrics  (METRIC). 

At  this  point,  the  code  is  fully  initialized  for  commencement  of  the  flow  field  sol- 
ution process,  which  is  conducted  via  an  iterative  loop.  At  each  iteration,  time  is  first 
marched  forward  one  time  step,  followed  by  computation  of  the  time  dependent  values 
of  angular  velocity,  angle  of  attack  and  step  change  of  angle  of  attack,  according  to  the 
relations: 

co  =  2kM00sm{2kM00[)  (2.20) 

a,=  a0  -  a,  cos(2A;Afoo0 
a,_,  =  a0  -  a,  cos[2AA/00(/  -  di)~\ 

da.  =  a,'  —  a,., 

The  grid  is  then  rotated,  in  the  physical  plane,  by  applying  the  the  following  relations 
at  each  grid  point. 
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x  =  x  cos(t/a)  —y  sin(t/a)  (2.21) 

y  =>•  cos((/a)  —  x  s'm(tisavepha) 

Following  recomputation  of  the  metrics,  the  solution  is  computed  by  subroutine  SLPS, 
the  ADI-algorithm,  which  calls  DISS1P,  for  computation  of  the  explicit  dissipation 
terms,  STRESS  and  RESI,  for  computation  of  the  inviscid  and  viscous  terms,  respec- 
tively, AMAT1  and  MATRIX!,  for  computation  of  the  Jacobian  and  its  inverse  in  the 
C  -direction  and  AMAT2  and  MATR1X2,  for  computation  of  the  Jacobian  and  its  in- 
verse in  the  r\  -direction.  Subroutine  STRESS  also  calls  subroutine  EDDY,  the  turbu- 
lence model,  for  computation  of  the  viscosity  coefficient.  WALLBC  enforces  the  wall 
boundary  conditions  and  finally,  solution  files,  as  discribed  later,  are  generated. 

Prior  to  resumption  of  the  loop,  performance  coefficients  are  generated  by  subrou- 
tines LOAD  and  CPPLOT.   Surface  pressure  coefficients  are  obtained  from  the  relation: 

CP  = Pb'P™        2  (2-22) 

l/2p00(A/00«00)2 

where  pb  is  the  surface  pressure  and  px  is  the  free  stream  pressure.  Skin  friction  coeffi- 
cients are  a  function  of  the  wall  shear  stress  (tJ  according  to  the  relations: 

Q= - r  (2.23) 

W=Uy-  Vx  £  J(^Xr  +  v^) 

The  aerodynamic  loads  of  lift,  drag  and  moment  about  the  quarter-chord,  are  computed 
by  the  following  relations. 

Q  =  Cn  cos(a)  —  C,  sin(a)  (2.24) 

Cd  =  C„  sin(a)  —  Ct  cos(a) 


C    = 


Q[<>  -ycu)dy  +  (x  -  xci*)dx^ 
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C„  and  C,  are  the  normal  and  tangential  forces  obtained  from  summing  surface  pressure 
and  skin  friction  forces  according  to  the  following. 


Cn 


Cpcix  + 


qdy^jc  {1.2: 


Q-[-  \cpdy+  fc/fr]/c 


G.     CODE  MODIFICATIONS 

The  following  modifications  were  made  to  the  code,  resident  on  the  FML  VAX, 
during  the  course  of  this  study. 

1.  Solution  output  commands  are  placed  within  the  loop  in  order  to  provide  interval 
output. 

2.  Read  commands  are  placed  at  the  beginning  of  the  main  program  in  order  to  allow 
storage  of  all  solutions  (previous  restarts  and  the  current  run)  in  a  sinde  combined 
file. 

3.  The  airfoil  section  is  rotated  to  the  initial  angle  of  attack,  vice  rotation  of  the  free- 
stream  direction. 

4.  Restarts  may  be  made  from  Plot3D  (output)  format  as  well  as  vector  format. 

5.  Additional  Plot3D  files  (surface  pressure  and  skin  friction  line  plots)  are  output 
from  subroutine  CPPLOT. 

6.  Grid  dimensions  are  made  true  variables. 

7.  User  inputs  are  reformatted  for  ease  of  use. 
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III.     DYNAMIC  GRAPHICS  GENERATION 

A.  INTERACTIVE  COMPUTER  GRAPHICS 

Due  to  the  complexities  and  time  dependent  phenomenon  associated  with  unsteady 
flows,  attainment  of  acceptable  clarity  in  flow  field  solutions  requires  complete  descrip- 
tive information  across  fine  grids,  at  a  large  number  of  minute  time  intervals.  Thus, 
analysis  and  visualization  of  data  generated  by  codes  such  as  Sankar's,  is  difficult  due 
to  both  the  volume  of  information  provided  and  its  temporal  nature.  Requirements 
identified  for  effective  How  field  visualization  include  maximization  of  the  bandwidth  of 
information  transfer,  such  that  it  closely  matches  the  capabilities  of  the  human  eye, 
maximization  of  the  quality  of  graphical  displays,  by  enhancing  key  features  and  sup- 
pressing others,  and  maximization  of  the  controllability  of  information  [Ref.  9].  Addi- 
tional advantages  in  information  transfer  can  be  achieved  through  the  use  of  redundant 
coding  and  structured  displays  [Ref.  10].  In  response  to  these  requirements,  interactive 
computer  graphics  workstations  have  been  evolved  to  complement  the  super-computer. 
Workstation  capabilities,  in  terms  of  geometrical  transformation  and  screen  update  dis- 
patch have  been  utilized  by  programmers  to  produce  effective  representations  of  flow 
field  motion,  through  synchronization  of  coordinated  data  sets.  Solution  clarity  is  en- 
hanced by  the  high  degree  of  spatial  resolution  and  large  range  of  colors  afforded  by  the 
workstation  displays.  The  interactive  capabilities  which  currently  exist  for  semi-complex 
flows  serve  to  improve  visual  cues  for  display  of  three-dimensional  data  sets  and  allow 
immediate  access  to  regions  of  interest. 

B.  HARDWARE 

The  NASA-ARC  Fluid  Mechanics  Laboratory  is  currently  equipped  with  an 
IRIS-3000  series  workstation  configured  with  4Mbytes  of  display  list  memory  and 
equipped  with  z-clipping  and  z-buffering  hardware  and  a  multiple  key  mouse.  The  IRIS 
display  processing  unit's  refresh  memory  is  arranged  as  a  two-dimensional  array  (768  x 
1024),  with  row-column  positions  matching  display  pixel  x-y  coordinates.  24  bits  are 
reserved  for  each  pixel  at  eight  bits  (256  intensities)  per  color,  resulting  in  a  range  of 
16,777,216  possible  colors.  When  displaying  dynamic  graphics,  the  refresh  memory  is 
divided,  in  order  to  meet  user  demands,  since  flyback  time  (about  1.3  n  sec)  is  too  short 
for  complete  picture  update.  This  double  buffering  mode  utilizes  half  the  memory  for 
refreshing  the  screen  displays,  while  the  other  half  is  dynamically  updated.   This  ensures 
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smooth  motion  on  the  display,  but  divides  the  number  of  bitplanes  per  pixel,  in  half. 
Thus,  the  number  of  colors  available  drops  to  4096.  In  order  to  avoid  a  similar  decrease 
in  the  range  of  colors  available,  color  maps  are  utilized.  In  this  case,  pixel  values  in  the 
refresh  memory  are  routed  to  an  index  or  color  look-up  table  (with  each  entry  composed 
of  24  predefined  bits,  which  define  the  color),  vice  direct  routing  to  the  intensity  digital- 
to-analog  converter.  The  IRIS  transformation  rate,  via  a  single,  composed  transforma- 
tion matrix  in  homogeneous  coordinates,  is  80,000  coordinates  per  second,  with  points 
and  lines  generated  at  3  million  pixels  per  second,  or  40,000  inches  per  second.  Polygons 
can  be  filled  (fiat  shading)  at  the  rate  of  44  million  pixels  per  second,  or  70,000  square 
inches  per  second.  These  capabilities  allow  generation  of  most  displays  at  the  rate  of 
approximately  10  screens  per  second,  which,  under  ordinary  conditions,  is  greater  than 
the  fusion  frequency  of  the  human  eye.  Interactive  transformations  are  ordinarily  ac- 
complished via  the  mouse. 

C.     PLOT3D  SOFTWARE 

Solutions  provided  by  the  Sankar  code  consist  of  multiple  binary  files,  specifically 
formatted  for  input  into  the  Plot3D  graphics  software,  authored  by  Pieter  Buning  for 
NASA.  Two  filetypes  exist,  which  are  labeled  as  XYZ  or  grid  files  and  Q  or  flow 
quantity  files.  Both  types  are  three-dimensional  matrices,  with  placement  of  data  within 
the  matrices  corresponding  to  the  row-column  addresses  of  mesh  points.  XYZ-files  have 
a  depth  of  two  in  the  two-dimensional  case  and  provide  cartesian  coordinate  definitions 
of  each  mesh  point.  Q-files  have  a  depth  of  four  in  the  two-dimensional  case  and  pro- 
vide computed  flow  quantities  at  each  mesh  point,  corresponding  to  the  q-vector  of 
equation  2.1.  Headers  for  each  file  list  free-stream  Mach,  Reynolds  number,  angle  of 
attack  and  time. 

Additional  flow  field  quantities  at  each  mesh  point  are  computed  according  to  the 
following  relations.  [Ref.  11] 

Definitions: 

y=  1.4 

R  =  1.0  (Gas  constant) 

M  =  Vic 


y£p 
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Based  on  these  values,  the  Plot3D  software  will  produce  plots  of  scalar  functions 
(density,  pressure,  temperature,  enthalpy,  energy,  velocity,  entropy,  momentum  and 
shocks)  suitable  for  contour  plots,  vector  functions  (velocity,  vorticity,  momentum  and 
pressure  gradient),  particle  trace  functions  (particle  traces  and  vortex  lines),  shock  wave 
locations  and  grid  functions  (computational  meshes  and  walls). 

Particle  traces  are  generated  using  trilinear  interpolation  of  values  of  the  vector 
function  inside  a  computational  cell  and  second-order  Runge-Kutta  steps  to  advance  the 
particle  in  space.  Five  steps  are  required  for  each  cell.  The  algorithm  for  computing 
shocks  computes  the  Mach  number  component  in  the  direction  of  the  local  pressure 
gradient.  Locations  where  this  value  decreases  through  1.0  is  plotted  as  a  shock.  Two- 
dimensional  stream  functions  are  calculated  by  integrating  the  mass  flow  across  a  coor- 
dinate line. 
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Output  from  Plot3D  is  available  in  a  variety  of  formats  and  is  a  function  of  user- 
defined  attributes.  For  dynamic  plots,  output  is  in  the  form  of  graphics  files,  formatted 
for  further  manipulation  by  animation  software.  Device  independent  plotting  (DIP)  is 
enabled  through  use  of  the  ARCGraph  binary  file  format  which  allows  data  to  be  ma- 
nipulated by  a  collection  of  function  libraries  and  utilities  developed  by  the  Advanced 
Computer  Research  Center  at  NASA- ARC.  Parameter  data,  along  with  graphics  prim- 
itives (device  independent  op-codes)  are  contained  in  these  files,  allowing  data  manipu- 
lation by  the  Cray,  IRIS  and  attached  VAXes.  Plot3D  generation  of  graphics  files 
requires  that  several  attributes  be  defined.  In  order  to  reduce  user-tasking,  Plot3D  will 
initially  search  for  a  filename-specific  initialization  communications  file,  which  contains 
all  required  inputs.  Input  files  must  be  read  singularly,  which  makes  necessary  the  sep- 
aration of  those  combined  plotting  files  received  from  the  Cray.  Output  graphics  files, 
however,  are  once  again  stored  in  combined  files  in  order  to  speed  further  processing. 
User-defined  attributes  include  axis  scales  and  extreme  values,  function,  number  of 
contours,  line  colors  and  types  and  wall  definitions. 

Subroutine  CPPLOT,  within  the  Sankar  code,  provides  grid  (XYZ)  files  which  con- 
tain plotting  information  for  surface  pressure  coefficient  and  skin  friction  coefficient  line 
plots,  corresponding  to  each  output  flow  field  solution  file.  These  line  plots  are  scaled 
and  translated  after  separation  from  the  combined  files  and  plotted  utilizing  the  grid 
function.  Likewise,  an  angle  of  attack  pointer  plot,  utilizing  a  sine  wave  format  (gener- 
ated external  to  the  code),  is  also  provided.  Q-files  containing  dummy  variables  (not 
utilized  for  grid  plot  generation)  are  also  provided  to  complement  the  line  plot 
XYZ-files,  as  required  by  Plot3D. 

D.     GRAPHICS  ANIMATION  SYSTEM  SOFTWARE 

The  final  step  in  the  visualization  process  is  the  animation  or  synchronization  of 
coordinated  graphics  datasets.  The  Graphics  Animation  System  (GAS)  is  a  software 
package  developed  by  Sterling  Software  under  contract  to  NASA  and  provided  to  edu- 
cational institutions  throughout  the  country.  It  is  device  specific,  running  only  on  IRIS 
workstations,  and  is  written  in  the  C  programming  language  under  the  UNIX  operating 
system.  Graphics  files,  in  ARCGraph  format,  are  read  and  stored  in  the  IRIS'  display 
list  memory  as  objects  by  the  menu-driven  program  and  assembled  according  to  se- 
quence file  instructions  (stored  transformation  matrices).  In  order  to  smooth  motion 
associated  with  geometric  transformations,  up  to  30  linearly  interpolated  matrices  are 
provided  by  the  program,  between  user-identified  keyframes.    Titles,  legends  and  object 
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attributes  are  also  available.  Ultimately.  How  Held  representations  may  be  transferred 
to  video  tape  or  16mm  film.  Interactive  viewing  (ordinarily  available  in  real-time)  is 
controlled  via  a  mouse,  within  the  menu's  "view  data"  window. 

A  complete  description  of  GAS  and  associated  video  hardware  may  be  found  in 
Reference  12. 
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IV.     RESULTS  AND  DISCUSSION 

A.  PROCEDURAL  CONSIDERATIONS 

Integration  of  the  F.ML  IRIS  with  related  CFD  software  and  the  Cray  X-MP/48 
was  tested  by  plotting  the  flow  field  about  an  airfoil  experiencing  deep  dynamic  stall,  as 
provided  by  the  Sankar  code.  The  advantages  of  plotting  dynamic  stall  are  two-fold. 
First,  it  is  an  extremely  complex  and  time-dependent  phenomenon  which  is  difficult  to 
visualize  from  discrete  data  sets.  Thus,  the  accuracy  of  solutions  for  this  type  flow  field 
can  only  be  determined  by  consideration  of  the  complete  flow  field  in  both  time  and 
space.  Dynamic  graphics  are  therefore  ideally  suited  for  its  study.  Secondly,  the  FML 
is  currently  heavily  involved  in  studies  of  dynamic  stall,  which  include  verification  of  the 
Sankar  code.    Cray  output  is  thus  used  to  full  advantage. 

In  order  to  provide  smooth  animations,  an  extremely  large  number  of  plotting  sets, 
provided  at  equal  time  intervals,  must  be  generated  by  the  code  during  the  oscillatory 
cycle.  Since  access  and  transfer  of  this  data  is  required  several  times  during  the  graphics 
generation  process,  it  was  found  that  storage  of  datasets  in  combined  files  provided  the 
greatest  efficiency,  especially  in  transfers  between  the  Cray  and  IRIS.  Software  resident 
on  the  IRIS  is  then  employed  to  separate  the  combined  file  into  individual  datasets. 
(Embedded  in  these  programs  are  schemes  to  scale  data,  as  necessary  to  control  their 
ultimate  on-screen  positions,  allowing,  for  example,  placement  of  line  plots  in  a  specific 
corner  of  the  display  or  plotting,  for  comparison,  two  flow  fields  side-by-side.)  Individ- 
ual datasets  are  automatically  assigned  names  which  correspond  to  those  contained  in 
Plot3D  initiation  files,  also  resident  on  the  IRIS.  These  files  currently  contain  com- 
mands for  the  processing  of  50  datasets.  Thus,  the  combined  file  separation  and 
ARCgraph  file  production  (Plot3X)  process  must  be  completed  at  intervals.  The  use  of 
repetitive  file  names  during  this  process,  however,  increases  automation  to  such  an  ex- 
tent that  only  a  limited  number  of  user  inputs  are  required  to  complete  the  process. 
ARCgraph  files  as  provided  by  Plot3D  are  again  in  combined  file  format,  allowing  easy 
(a  single  user  command)  input  into  animation  software. 

B.  THE  DYNAMIC  STALL  PROCESS 

A  prerequisite  to  review  of  any  code-provided  flow  field  solution  is  consideration  of 
available  empirical  information.  The  following  dynamic  stall  characteristics  have  been 
observed  for  both  harmonically  oscillating  airfoils  [Ref.   13]  and  airfoils  undergoing 
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monotonic  angle  of  attack  increases  [Refs.  14,  15].  First,  the  airfoil  stalls  at  an  angle  of 
attack  which  exceeds  the  static  stall  angle.  Second,  the  corresponding  normal  forces  and 
pitching  moments  exceed  those  attainable  in  the  static  case  and,  finally,  a  rapid  change 
in  pitching  moment  magnitude,  termed  as  "moment  stall"  by  Harris  and  Pruyn  [Ref.  16] 
occurs  several  degrees  in  azimuth  prior  to  the  normal  force  decrease,  or  "lift  stall",  vice 
simultaneous  occurence  as  in  the  static  case.  The  process  which  results  in  these  char- 
acteristics can  be  broken  down  into  a  series  of  chronological  events,  as  depicted  in  Fig- 
ure 2  from  Reference  13.  While  the  specific  characteristics  for  a  given  airfoil  are  a 
function  of  free-stream  Mach  number,  Reynolds  number,  reduced  frequency,  mean  angle 
of  attack  and  oscillation  amplitude,  empirical  data  obtained  from  a  NACA-0012  airfoil 
at  k  =  .15,  Re  =  2.5xl06  and  a  =  151  —  10'  sin(a)f)  is  provided  in  the  following  discussion 
to  illustrate  the  dynamic  stall  process. 

At  point  (a)  of  Figure  2,  the  static  stall  angle  is  exceeded  without  any  detectable 
chanse  in  flow  over  the  airfoil.  The  boundarv  laver  remains  thin  with  no  evidence  of 
flow  reversal  at  the  surface.  At  point  (b),  (  a  =  19  —  20')  flow  reversal  appears  at  the 
surface.  Over  the  majority  of  the  airfoil,  the  boundary  layer  remains  thin  and  attached. 
At  the  rear  portion  of  the  airfoil,  however,  the  boundary  layer  thickens  as  the  flow 
gradually  decreases  to  zero  velocity.  At  point  (c),  large  eddies  appear  in  the  boundary 
layer.  By  point  (d),  flow  reversal  has  spread  over  much  of  the  chord,  from  the  trailing 
edge  to  x/c  =  0.3.  With  as  much  as  50%  of  the  airfoil  experiencing  flow  reversal,  no 
changes  in  C„  or  Cm  are  detected.  At  point  (e),  (  a  =  23.4°  )  a  vortex  forms.  The 
boundary  layer  at  the  front  of  the  airfoil  abruptly  breaks  down  (simultaneously  from  x  c 
=  0.0  to  0.3)  and  the  vortex  moves  downstream  at  approximately  35  -  45%  of  the  free- 
stream  velocity.  At  point  (0  the  lift-curve  slope  increases  beyond  the  Ina.  limit  of 
quasi-static  flows.  By  point  (g),  the  pressure  distribution  is  altered  sufficiently  to 
produce  a  noticeable  divergence  in  the  pitching  moment  (moment  stall)  but  is  accom- 
panied by  a  continued  increase  in  lift.  Maximum  lift  occurs  at  approximately  mid-chord, 
followed  by  a  sharp  decrease  (point  (h)),  which  identifies  lift  stall.  (Boundary  layer  sep- 
aration has  occurred  to  such  an  extent  that  continued  increases  in  angle  of  attack  will 
not  result  in  continued  increases  in  lift.)  At  point  (i)  (  a  =  24.95  °  ),  maximum  negative 
moment  occurs.  The  vortex  moves  off  the  trailing  edge  at  a  =  24.8  °  (downstroke)  and 
C„  and  Cm  move  toward  static  levels.  At  point  (j),  the  airfoil  is  fully  stalled.  At  point  (k), 
boundary  layer  flow  begins  to  reattach,  progressing  from  the  leading  edge  at  approxi- 
mately 25  -  35%  free-stream  velocity  and  is  complete  by  a  =  7  °  (downstroke).  At  point 
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(a)   STATIC  STALL  ANGLE  EXCEEDED 


(b)   FIRST  APPEARANCE  OF  FLOW 
REVERSAL  ON  SURFACE 


(c)    LARGE  EDDIES  APPEAR  IN 
BOUNDARY  LAYER 


10  15  20 

INCIDENCE,  a,  dag 


(d)    FLOW  REVERSAL  SPREADS  OVER 
MUCH  OF  AIRFOIL  CHORD 


(e)   VORTEX  FORMS  NEAR 
LEADING  EDGE 


(f)    LIFT  SLOPE  INCREASES 


(g)   MOMENT  STALL  OCCURS 


(h)    LIFT  STALL  BEGINS 


(i)   MAXIMUM  NEGATIVE  MOMENT 


(j)    FULL  STALL 


(k)    BOUNDARY  LAYER  REATTACHES 
FRONT  TO  REAR 


(I)    RETURN  TO  UNSTALLED  VALUES 


Figure  2.      The  Dynamic  Stnll  Process  (from  Carr,  Ref.  13) 
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(1),  after  some  hysteresis  (  a  =  6  °  ,  upstroke)  the  separated  region  has  fully  closed  and 
unstalled  values  oi~  C„  and  Cm  are  reestablished. 

C.     PROCEDURAL  VERIFICATION 

The  case  chosen  for  procedural  verification  closely  corresponds  to  the  Carr  data  of 
the  previous  section  and  exactly  matches  conditions  previously  run  when  the  code  was 
initially  vectorized  for  use  on  the  ARC  Cray  [Ref.  17].  There  are  several  advantages 
associated  with  this  duplication.  First,  it  allowed  storage  of  a  complete  record  of  this 
case  (data  saved  at  300  intervals,  vice  12)  for  future  access.  Second,  it  provided  a  simple 
means  by  which  to  determine  if  code  modifications  were  correctly  incorporated.  Finally, 
it  provided  a  baseline,  in  conjunction  with  empirical  data,  against  which  future  sensitiv- 
ity studies  could  be  compared. 

Inputs  for  this  case  were  as  follows: 

Re  =  3.45xl06 

AC  =  -2S3 

6t  =  .005  sec. 

>1  -nun  =  .00005 

co  =  .151 

Artificial  Viscosity  =  5 

Grid  size:  157  x  40 

Oscillatory  Cycle:  a  =  15°  -  10°  cos(cor) 

The  original  Baldwin-Lomax  turbulence  model  was  used  throughout  the  oscillatory 
cycle.  Figures  3  through  7  are  coefficient  of  pressure  contour  plots  of  the  predicted  flow 
field.  (Dynamic  plots  are  similar  to  these  Plot3D-provided  plots,  but  do  not  contain 
explicit  quantitative  information.)  As  previously  noted  [Ref.  17,  p.  48],  under  these 
conditions,  moment  coefficients  are  consistently  underpredicted  as  compared  to  exper- 
imental data.  Drag  and  lift  coefficients  closely  match  experimental  data  below  approxi- 
mately 15  and  18  degrees  angle  of  attack,  respectively.  Early  prediction  of  flow 
separation,  however,  forces  inaccuracies  in  quantitative  solutions  beyond  these  values 
and  underprediction  of  the  maximum  lift  coefficient  obtainable.  These  findings  highlight 
the  rationale  behind  the  modified  turbulence  model.  Information  provided  by  the  dy- 
namic plots,  shows  general  qualitative  agreement  with  experimental  data,  during  the 
upstroke.  A  vortex  forms  near  the  leading  edge  of  the  airfoil  section  and  progresses 
down  its  upper  surface,  gradually  increasing  in  size.  As  the  downstroke  begins,  however, 
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Figure  3.       Pressure  Contours,  a  =   15°  -  10°  (.15 It),  a  =  5.00°,  5.24  °,  5.92°,  7.02°, 
8.47°,  10.23°  ,  (top  to  bottom,  left  to  right),  Mx  =  .283,  Re  =  3.45  x  10\ 
r]mm  =  .00005. 
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Figure  4.  Pressure  Contours,  a  =  15°  -  IU°  (.151t),  a  =  12.20°,  14.30°,  16.43°, 
18.50°,  20.40°,  22.07°  ,  (top  to  bottom,  left  to  right),  M,  ~  .283, 
Re  =  3.45  x  IU\  >/mjn  =  .00005. 
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Figure  5.  Pressure  Contours,  a  =  15°  -  10°  (.15lt),  a  =  23.41°,  24.36°,  24.89°, 
24.98°,  24.60°,  23.80°  ,  (top  to  bottom,  left  to  right),  A/..  =  .283, 
fo  =  3.45  x  10%  rjmin  =  .00005. 
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Figure  6. 


Pressure  Contours,  a  =  15°  -  10°  (.151t),  a  =  22.59°,  21.03°,  19.20°, 
17.19°,  15.07°,  12.9-4°  ,  (top  to  bottom,  left  to  right),  MM  =  .283, 
fo=  3.45  x  10*,  Vmm  =  . 00005. 
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Figure  7.       Pressure  Contours,  a  =   15°  -  10°  (.1511),  a  =   10.92°,  9.07°,  7.50°,  6.27°, 
5.43°,  5.00°  ,  (top  to  bottom,  left  to  right),  M„  =  .283,  K<z  =  3.45  x  10\ 
>/m,n  =  -00005. 
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rather  than  shedding  from  the  trailing  edge,  the  vortex  stagnates  at  approximately  the 
65%  chord  point  and  gradually  dissipates.  Line  plots  of  surface  pressure  coefficients 
more  accurately  portray  (as  compared  to  traditional  carpet  plots)  the  turbulent  effects 
associated  with  the  existence  of  the  vortex.  During  the  vortex  dissipation  phase,  these 
line  plots  show  a  slight  (gradually  decreasing)  pressure  rise  in  the  vicinity  of  the  vortex. 
The  subtle  differences  between  these  plots  and  those  which  would  have  been  obtained 
had  the  vortex  actually  shed,  are  lost  when  this  information  is  integrated  for  lift,  drag 
and,  to  a  lesser  extent,  moment  coefficient  data.  Thus,  quantitative  data  may  not  be 
eood  indicators  of  flow  field  solution  accuracy. 

D.     SENSITIVITY  ANALYSIS 

The  following  code  parameters  which  should  effect  dissipation  of  the  vortex  have 
been  identified. 

1.  Artificial  viscosity  magnitude. 

2.  Turbulence  model  parameters. 

3.  Grid  resolution. 

A  limited  sensitivity  analysis  into  the  effects  of  grid  resolution  changes  in  the  trailing 
edge  region  was  conducted  in  order  to  provide  direction  for  follow-on  studies.  During 
the  initial  phase  of  this  analysis,  only  qualitative  results  will  be  considered,  for  those 
reasons  mentioned  above.  Dynamic  graphics,  therefore,  are  ideally  suited  for  this  task, 
especially  when  data  scaling  capabilities  are  utilized  in  order  to  plot  comparative  data- 
sets,  side  by  side. 

The  grid  is  relatively  coarse  in  the  region  in  which  the  vortex  becomes  stationary  and 
dissipates,  as  compared  to  those  regions  in  which  it  behaves  as  expected.  In  the  first 
sensitivity  case,  therefore,  grid  dimensions  were  held  constant  and  the  distance  of  the 
first  constant-  r\  line  from  the  airfoil  surface  was  increased  to  .001  ,  resulting  in  a  finer 
grid  in  the  trailing  edge  region.  Figure  S  shows  the  effects  of  this  change  on  the  vortex. 
Within  this  grid,  the  vortex  once  again  becomes  stationary,  but  effectively  splits,  such 
that  a  secondary  vortex  is  shed  while  the  original  vortex  dissipates.  While  this  does  not 
correspond  to  any  process  encountered  in  actual  flows,  it  does  indicate  that  final  sol- 
utions are  sensitive  to  this  parameter. 

Changes  in  r\  -spacing  have  two  effects  on  the  flow  field  solution.  Grid  spacing  in 
the  trailing  edge  region  becomes  more  fine  at  the  expense  of  grid  spacing  in  the  boundary 
layer,  which  becomes  more  coarse.    The  extent  of  influence  of  the  boundary  layer  sol- 
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Figure  8.  Pressure  Contours,  a  =  15°  -  10°  (.15 It),  a  =  24.89°,  24.85°,  23.80°, 
21.85°,  19.20°,  16.13°  ,  (top  to  bottom,  left  to  right),  M^  =  .283, 
Re  =  3.45  x  10*,  rjmin  =  .001. 
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ution  on  the  vortex,  or  far  field  solution,  is  not  yet  clear.  It  is  known,  however,  that 
boundary  layer  solutions  are  extremely  sensitive  to  ;/  -spacing  parameter  changes,  re- 
quiring an  initial  value  of  .00005  for  accurate  solutions  [Kef.  IS],  Further  study,  there- 
fore, is  required  in  order  to  determine  whether  vortex  behavior  changes  observed  in  the 
previous  case  were  in  response  to  the  altered  boundary  layer  solution  (which  would  re- 
commend investigation  of  turbulence  model  parameters)  or  the  altered  grid.  In  order  to 
facilitate  this  study,  the  existing  grid  would  require  enlargement  (from  161  x  4h.  Un- 
fortunately, this  results  in  such  a  large  memory  allocation  from  the  Cray,  that  it  is  pro- 
hibitively inefficient.  Replacement  of  the  entire  grid  generation  process  is  therefore 
currently  under  consideration  for  this  code. 

E.     CONCLUSIONS 

The  current  study  has  resulted  in  completion  of  the  following  items. 

1.  Full  procedural  documentation  (Appendix  A)  has  been  developed  for  efTieient, 
code-independent,  two-dimensional,  dynanuc  graphics  production  on  FML  and 
NTS  IRIS  workstations.  Data  must  be  in  Plot3D  format  and  may  be  provided  by 
either  code  or  experiment.  A  method  for  simultaneous  viewing  of  multiple  datasets 
has  been  incorporated. 

2.  The  Sankar  Navier-Stokes  solver  (Appendix  B)  has  been  modified  io  allow  dynamic 
graphical  presentation  of  its  flow  field  solutions. 

3.  Dynamic  graphics  applications  in  the  study  of  complex  flows  were  illustrated  by 
means  of  an  introductory  sensitivity  analysis,  using  the  Sankar  code.  This  limited 
analysis  identified  either  grid  resolution  in  the  trailing  edge  region  or  the  boundary 
layer  solution  as  parameters  to  which  vortex  behavior  is  strongly  dependent.  (This 
suggests  that  incorporation  of  a  more  sophisticated  grid  generation  scheme  would 
be  worthwhile.) 

In  order  to  derive  full  benefit  from  application  of  this  technology,  existing  and  future 
empirical  data  should  be  stored  in  plotting  format,  allowing  the  availability  of  visual  re- 
cords of  multiple  flow  functions.  Also,  modification  of  additional  Navier-Stokes  solvers 
for  dynamic  output  will  allow  effective  code-to-code  and  code-to-experiment  compar- 
isons, all  in  an  effort  to  eventually  develop  an  accurate  flow  field  predictor. 
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APPENDIX  A.     PROCEDURES  FOR  GENERATING  DYNAMIC 

GRAPHICS  AT  FML 

A.     THE  ANIMATION  PROCESS 

Creation  of  dynamic  graphics  for  flow  field  visualization  at  the  NASA-Ames  Re- 
search Center  Fluid  Dynamics  Laboratory  involves  five  basic  processes: 

1.  Navier-Stokes  solver  is  subn  itted  to  the  Cray  X- MP/48  and  dynamic  output  saved 
in  the  Central  Storage  Facility  (CSF)  in  combined  files. 

2.  Combined  files  are  converted  to  IRIS  binary  and  transferred  to  an  FMLIRISl  ac- 
count. 

3.  Combined  files  are  separated  into  individual  plotting  files  on  the  IRIS  and  con- 
verted to  combined  graphics  (.gra)  files  by  Plot3D. 

4.  Graphics  files  are  fed  to  the  Graphics  Animation  System  software  (GAS)  and  dy- 
namic graphics  are  plotted. 

5.  Resultant  files  and  plots  are  transferred  to  storage  or  video  tape. 

Prior  to  beginning  the  animation  process,  the  code  must  be  modified  to  provide  dy- 
namic output  in  Plot3D  format.   This  involves: 

1.  Placement  of  solution  output  commands  within  the  iterative  loop  along  with  some 
interval  flag. 

2.  If  restarts  are  to  be  used,  placement  of  read  statements  at  the  beginning  of  the 
program  to  allow  storage  of  all  provided  solutions  in  a  single  combined  file. 

3.  Storage  of  any  additional  required  plotting  information  (line  plots)  in  Plot3D  XYZ 
files  for  eventual  display  utilizing  Plot3D's  grid  function. 

Full  documentation  of  Plot3D  formats  and  functions  may  be  found  in  Sterling  Software 
technical  note  15,  A  Guide  to  Generating  Movies  using  PlotlD  and  GAS,  by  Greg  Howe. 
The  procedures  outlined  below  refer  specifically  to  a  version  of  the  Sankar  Navier-Stokes 
solver,  modified  for  dynamic  output.  They  can  easily,  however,  be  generalized  to  apply 
to  any  code  (or  experiment)  from  which  dynamic  output  is  desired.  While  many  meth- 
ods of  completing  these  steps  may  exist,  those  outlined  below  have  proven  to  be  the 
most  time  efficient. 

1.     Vax  Submitted  Cray  Jobs 

A  full  dynamic  cycle,  utilizing  a  161  x  41  grid,  requires  over  5000  seconds  of 
Cray  time,  making  job  chaining  necessary  to  obtain  a  reasonable  priority.  With  900 
second  jobs,  24  to  48  hours  are  usually  sufficient  for  the  full  computation  to  be  com- 


pleted.  Since  significant  differences  exist  between  the  first  JCL  from  which  Plot3D  tiles 
are  saved  and  subsesquent  restart  JCL's,  two  JCL's  are  maintained  on  identical  copies 
of  Sankar's  code. 

INITIAL. JCL.  Accesses  CSF  or  Cray  tape  (via  STAG  EX  commands)  for  initial 
solution,  which  may  be  in  either  Plot3D  or  matrix  format.  Combined  Plot3D  files  are 
initialized  on  CSF  and  assigned  PDN's  and  IDs,  Tape  05  input  data  (appended  to  end 
of  code)  is  assigned  in  accordance  with  definitions  included  in  MAIN  portion  of  the 
code.  (If  restart  is  from  dynamic  data,  AOA  will  not  be  a  whole  number.  TSTART 
must  be  adjusted  accordingly,  in  order  to  have  accurate  time  AOA  matches.) 

RESTART. JCL.  Tape  05  inputs  must  identically  match  those  from 
INITIALJCL,  with  the  exception  of  TSTART,CSTP,CPLT  and  possibly  FORMAT. 
INITIAL. JCL-defined  PDN  ID  solutions  are  accessed  (read,  stored  and  appended  with 
additional  solutions),  saved  (with  identical  PDN  IDs)  and  scrubbed  (oldest  versions 
deleted).  The  initial  restart  solution  is  accessed  with  library  routine  GETP3D.  written 
by  Greg  Howe,  which  will  read  a  combined  file,  skip  a  specified  number  of  datasets 
(DSSKIP)  and  access  the  next  dataset  (or  datasets,  if  NUMDS  is  not  equal  to  one,  the 
default).    Thus,  if 

DSSKIP  +  CPLT-  1, 

the  correct  solution  will  be  provided  for  restarts.  Subsequent  restarts  illustrate  the  ad- 
vantage of  combined  file  name  (PDN  ID)  repetition,  as  they  simply  require  the  following 
RESTART. JCL  and  Tape  05  updates: 

1.  DSSKIP 

2.  Job  chain  commands. 

3.  CSTP 

4.  CPLT 

Once  the  optimum  number  of  plotting  sets  for  a  cycle  is  determined,  plotting 
files  for  an  AOA  pointer  on  IRIS  displays  can  be  generated  using  Plot3D.  (500  sets 
maximum,  or  modify  dimension  statements.) 

Line  Plots.  Any  pointer  or  line  plot  information  may  be  generated,  either  within 
the  code  or  externally,  and  stored  in  XYZ  files  for  further  processing.  Utilization  of 
Plot3D  function  0  and  proper  descriptors  of  walls  will  produce  the  desired  plots. 
SINE. TXT  is  a  program  which  generates  plots  for  a  dynamic  AOA  pointer  in  sine  wave 
format.    In  this  program,  Tape  05  includes  values  for  pointer  scaling  (axis  length  with 
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no  scaling  is  In  )  and  positioning  on  the  IRIS  screen  and  initial  pointer  location.    For 
example,  if  the  cycle  starts  from  the  minimum  AOA, 

PNLOCI  =  .75(NPLOTS). 

Hardcopy  printouts  may  also  be  obtained  from  these  grid  files  by  use  of  pro- 
grams such  as  CPPLOT.TXT.  This  program  pro\ides  printouts  of  upper  and  lower 
surface  pressure  coefficient  and  skin  friction  values  at  x  c  locations,  including  Co  plot 
generation. 

2.  Path  to  the  IRIS 

SEND.JCL.  Utilizes  GETP3D  and  SENDWKS  to  convert  to  IRIS  binary  and 
transfer  complete  or  partial  combined  files,  or  individual  plotting  sets,  from  the  CRAY 
to  the  IRIS.  (Since  I  O's  between  the  Cray  and  IRIS  can  get  "clogged'  if  too  many 
transfers  are  requested,  the  individual  method  is  not  recommended.)  This  may  need  to 
be  completed  in  stages  to  avoid  exceeding  the  IRIS'  memury  capacity. 

3.  Creating  Graphics  Files 

Combined  files  on  the  IRIS  are  separated  into  individual  plotting  sets  by  the 
programs  GETSIN.F  (also  generates  the  dummy  Q  tile,  "qsin.iris"),  GETCP.F  (also 
scales  and  positions  Cp  or  skin  friction  on  the  IRIS  screen)  and  GETX.F  (for  XYZ  and 
Q  flow  field  datasets).  These  programs  require  some  initial  user  inputs  prior  to  running. 
Output  dataset  names  will  automatically  match  Plot3D  initialization  file  names.  Since 
the  Plot3x  run  for  each  dataset  type  (X,  P  and  SIN)  requires  specific  arguments,  initial- 
ization files  for  each  data  type  exist  separately  in  files  XINI.COM,  PINT.COM  and 
SINT.COM.  Since  Plot3X  searches  for  the  file  name  "PLOT3DINT.COM',  these  flies 
must  be  renamed  appropriately,  utilizing  the  "mv"  command.  Entering  the  command 
"plot3x"  will  then  generate  multiple  graphics  files  which  are  stored  in  a  combined  file 
matching  the  initial  Q  file  name  encountered  by  Plot3D,  (i.e.,  "qUOl.gra").  This  file  is 
then  stored  on  an  appropriate  subdirectory  or  fed  to  GAS. 

4.  Notes  on  GAS 

For  movies,  the  maximum  number  of  objects  per  sequence  is  fifty,  The 
"seq(n).seq"  files  are  50  object  far-field  solution  files  (seq(n).seq  =  objects  n(l-5G|), 
which  use  object  400  for  titles.  They  can  be  fed  to  GAS  (after  object  input)  via  the  AL'X 
10  window.  For  interactive  viewing,  the  VIEW  DATA  window  is  the  only  usable  op- 
tion. (No  sequence  files  are  required  in  view  data.)  Full  GAS  documentation  is  avail- 
able in  the  GAS  User's  Guide,  available  from  Stei  Una  Software. 
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5.     Storage 

Since  IRIS  memory  space  is  severely  limited  at  FML,  resultant  graphics  files 
should  be  stored  elsewhere,  if  not  in  immediate  use.  Downloading  to  1  4"  cassette  tape 
is  easily  accomplished  (also  a  recommended  backup),  as  is  transfer  to  CSF. 
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B.     SUPPORT  CODES 
1.      INITIAL.JCL 


JOB  ,  JN-FLYNAVY  ,  T-900  . 
ACCOUNT  ,  AO .  US-  .  UPW-  . 
•  . 
CFT.ON-A.OFF-S. 


ERIC  PAGENKOPF.X4269. 


•RESTART  COMMANDS: 

CSF .  ACCESS  .  CSFACCT- . CSFPSWD- . DN-OLDSLN ,  PDN-NSE509 .  I D-VALDES . 
ASS  IGN . DN-OLDSLN . A-FT07 . 


•.SAVE  DATA  FROM  PRESENT  RUN  IN  COMBINED  FILE: 
ASSIGN.DN-XYZ.    A-FT20. 
ASSIGN.DN=Q.      A-FT21. 
ASSIGN.DN-CPX.    A-FT50. 
ASSIGN. DN-CFX,    A-FT60. 


LDR , MAP-PART , SET-ZERO . 

•  . 

CSF . SAVE . CSFACCT- . CSFPSWD- . DN-XYZ , 

CSF . SAVE . CSFACCT- . CSFPSWD- . DN-0 , 

CSF . SAVE . CSFACCT- . CSFPSWD- . DW-CPX . 

CSF . SAVE . CSFACCT- , CSFPSWD- , DN-CFX . 


PDN-SNKR05X , 
PDN-SNKR05O, 
PDN-SNKR05P . 
PDN-SNKR05F. 


ID-FLYNAVY. 
ID-FLYNAVY. 
ID-FLYNAVY. 
ID-FLYNAVY. 


. ACCESS . DN=SENDVAX . PDN-SENDVAX . ID-STTRDM . OWN-RFTRDM . 

.SENDVAX. DN-XYZ. VDN-'RAL" J IANPS  password" :  : [ J  I ANPS . PAGAN. DATA]XT1 .DAT' . 
. SENDVAX.DN-Q, VDN='RAL"J IAH'S  password" : : [J IANPS. PAGAN. DATA]QT1 .DAT" . 
. ACCESS . DN-SENDWKS , PDN-SENDWKS . ID-STTRDM . OWN-RFTRDM . 

SENDWKS.DN=03X49.VDN=,FMLlRISr,j  i ao  password": : "/u/j  iaa/pagan/q.  iris"' .NOREC. 
.SENDWKS. DN-CFX. VDN-' FMLIRIS1" j i aa  password": : "/u/j i aa/pagan/cf 01 .dot"' .NOREC. 


JOB  CHAINING  COMMANDS: 

FETCH . DN-JOB2 . TEXT- ' [ PAGANjNSMULT . TXT ; 2 ' 

REWIND.DN-JOB2. 

SUBMIT.DN-JOB2. 


/EOF 
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REST  ART.  J CL 


JOB . JN-FLYNAVY . T-900 . 
ACCOUNT  .  AO ,  US» .  UPW- . 
•  . 
CFT.ON-A.OFF-S. 


ERIC  PAGENKOPF.X4269. 


•.RESTART  COMMANDS: 

•  .  CSF . ACCESS . CSFACCT- . CSFPSWD- . DN-OLDSLN . PDN-NSE509 
CSF . ACCESS . CSFACCT- . CSFPSWD- . DN-GETP3D ,  PDN-GETP3D . 


CSF , ACCESS . CSFACCT- . CSFPSWD= . DN=XYZ 1 . 

CSF, ACCESS. CSFACCT-. CSFPSWD- . DN-Q1 , 

CSF , ACCESS . CSFACCT- . CSFPSWD- . DN-CPX1 . 

CSF . ACCESS . CSFACCT- . CSFPSWD- . DN-CFX 1 , 

REWIND. DN-Q1 . 

GETP3D. I-Q1 .O-OLDSLN.DATATYPE-Q.DSSKIP-49 

REWIND. DN-OLDSLN. 

REWIND.DN-XYZ1 . 

REWIND.  DN-Q1  . 

REWIND.DN-CPX1 . 

REWIND.DN=CFX1  . 

ASSIGN.DN-XYZ1 .   A-FT90. 

ASSIGN.DN=01 ,     A-FT91. 

ASSIGN.DN-CPX1 .   A-FT92. 

ASSIGN,DN=CFX1 .   A-FT93. 

ASS  I GN . DN=OLDSLN . A-FT07 . 


PDN=SNKR05X. 
PDN=SNKR05Q. 
PDN=SNKR05P. 
PDN-SNKR05F. 


ID-VALDES. 

ID-RFRICC. 

ID-FLYNAVY. 

ID-FLYNAVY. 

ID-FLYNAVY. 

ID-FLYNAVY. 


..SAVE  DATA  FROM  PRESENT  RUN  IN  COMBINED  FILE:- 
ASSIGN.DN-XYZ.    A-FT20. 
ASSIGN.DN=Q.      A-FT21. 
ASSIGN.DN-CPX.    A-FT50. 
ASSIGN.DN-CFX.    A-FT60. 


LDR . MAP-PART , SET-ZERO . 


CSF, SAVE. CSFACCT- 
CSF. SAVE .CSFACCT- 
CSF. SAVE. CSFACCT- 
CSF.  SAVE. CSFACCT- 
CSF. SCRUB. CSFACCT 
CSF, SCRUB. CSFACCT 
CSF. SCRUB. CSFACCT 
CSF. SCRUB .CSFACCT 


PDN-SNKR05X. 
PDN-SNKR05Q. 
PDN-SNKR05P. 
PDN-SNKR05F. 


CSFPSWD-.DN-XYZ. 
CSFPSWD-. DN-Q. 
CSFPSWD-.DN-CPX. 
CSFPSWD-.DN-CFX. 
=, CSFPSWD-. PDN-SNKR05F, ID-FLYNAVY. 
.CSFPSWD-.PDN-SNKR05P. ID-FLYNAVY. 
,CSFPSWD=.PDN=SNKR05X, ID-FLYNAVY. 
. CSFPSWD- . PDN-SNKR05Q . I D-F  LYNAVY . 


ID-FLYNAVY. 
ID-FLYNAVY. 
ID-FLYNAVY. 
ID-FLYNAVY. 


. ACCESS . DN-SENDVAX . PDN-SENDVAX .ID-. OWN- . 

. SENDVAX.DN=XYZ.VDN-'RAL"JIANPS  possword" :  : [ J  I ANPS . PAGAN. DATA]XT1 .DAT" . 

.SENDVAX,DN=Q,VDN=-RAL"J1ANPS  pos sword": : [ J 1ANPS. PAGAN . DATA] QT1 .DAT" . 

. ACCESS . DN-SENDWKS . PDN-SENDWKS . I D= . OWN- . 

. SENDWKS.DN=Q3X49.VDN="FMLIRIS1" j >oa   password": : "/u/j  iao/pagan/q.  iris"" .NOREC, 

.  SENDWKS.DN=CFX,VDN=,FMLIRIS1" j ioa  password" : :"/u/]  i aa/pagon/cf 01 .dot'"  , NOREC. 


JOB  CHAINING  COMMANDS 


FETCH . DN-JOB3 . TEXT- ' [PAGAN] JOB3 .  TXT ' 

REWIND.DN-JOB3. 

SUBMIT.DN-JOB3. 


/EOF 
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3.   SEND.JCL 


JOB, JN-GETX.T-30.  ERIC  PAGENKOPF,  X4269 . 
ACCOUNT , AC= , US= . UPW- . 

•  . 

CSF . ACCESS . CSFACCT= , CSFPSWD- , DN=XYZ . PDN=SNKR38X . I D-FLYNAVY . 
CSF , ACCESS , CSFACCT= , CSFPSWD- . DN-Q . PDN-SNKR38Q . I D-FLYNAVY . 

•  . 

CSF .  ACCESS .  CSFACCT= ,  CSFPSWD= ,  DNNGETP3D ,  PDN-GETP3D  ,  I D-RFRICC . 
ACCESS . DN=SENDWKS . PDN=SENDWKS . I D= ,OWN= . 

•  . 

REWIND, DN=XYZ. 
REWIND, DN=Q. 

• 

GETP3D. I=XYZ.O=XOUT1 .DATATYPE-XYZ.DSSKIP-100 .NUMDS-150. 
GETP3D, I-O.0=OOUT1 ,DATATYPE=Q.DSSKIP=100. NUMDS=150. 

• 

SENDWKS.DN=XOUT1 . VDN= ' FML IRIS1 " j iaa  possword":  : "/u/j i aa/pogan/38x . iris"'  , NOREC. 

SENDWKS.DN=OOUT1 , VDN= ' FML IRI S1 " j ioa  password": : "/u/j  i aa/pagan/38q. iris"1 .NOREC. 

•  . 

• . FETCH. DN=JOB2,TEXT=" [PAGAN ] SEND. JCL;2' . 

• . REWIND. DN=JOB2. 

•.SUBMIT.DN=JOB2. 
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4.   SINE.TXT 


JOB, JN-FLYNAVY.T-30.   ERIC  PAGENKOPF, X4269 . 
ACCOUNT  ,  AO .  US-  .  UPW- . 

•  . 

ACCESS . DN-SENDWKS . PDN-SENDWKS .ID-. OWN- . 

CFT.ON-A.OFF-S. 

•  . 

ASSIGN. DN-X.   A-FT50. 
ASSIGN. DN-Q,   A-FT60. 

• 

LDR . MAP-PART , SET-ZERO . 

• 

SENDWKS.DN=0,VDN='FMLIRIS1" j  i  oo  PASSWORD":  :/u/j  i  oo/pogan/qa  i  n  .  iris'  .NOREC. 
SENDWKS.DN-X.VDN-'FMLIRISI" j  i oa  PASSWORD":  :/u/j  i aa/pagan/snk r05s '  .NOREC. 
CSF , SAVE . CSFACCT- . CSFPSWD- , DN=X . PDN-SNKR05S . I D-FLYNAVY . 


/EOF 


PROGRAM  POINTER 


THIS  PROGRAM  GENERATES  PL0T3D  FILES  WHICH  WILL  PRODUCE  AN  • 

AOA  POINTER  WHEN  FUNCTION  0  IS  SELECTED.  XY2  FILES  ARE  SAVED  • 

ON  CSF  FOR  FUTURE  ACCESS  VIA  PROGRAM  GETSIN.JCL.  A  SINGLE  Q  • 

FILE  OF  PROPER  DIMENSION  IS  SENT  DIRECTLY  TO  THE  IRIS  WITH  • 

FILENAME  QSIN.IRIS.  • 

• 

VARIABLES:  • 

NPLOTS:  TOTAL  NUMBER  OF  PLOTS  IN  ONE  CYCLE.  • 

PNLOCI:  INITIAL  INTEGER  LOCATION  OF  POINTER  • 

PNLOC:  INTEGER  LOCATION  OF  POINTER.  • 

AXIS  SCALING  FACTOR  • 

SINE  WAVE  SCALING  FACTOR  • 

X-POSITIONAL  FACTOR  (SCREEN  LOCATION)  • 

Y-POSITIONAL  FACTOR  (SCREEN  LOCATION)  • 

• 


•  XSCALE 

•  YSCALE 

•  XPOSIT 

•  YPOSIT 


INTEGER  PNLOC. GDIM. PNLOCI 

DIMENSION  SINX(3.500) .SINY(3,500) .01 (3 .500) .02(3 . 500) 
DIMENSION  03(3.500) .04(3.500) 
DATA  01  . 02.03.04/6000*1  ./ 
C 

C««»   INITIALIZE 
C 

READ  (5.1) 

READ  (5.2)  NPLOTS. PNLOCI .XSCALE. YSCALE. XPOSIT. YPOSIT 
PI=ATAN(1 . )»4. 
GDIM  -  3 
PNLOC  -  PNLOCI 
C 

DO  10  L-1 .NPLOTS 
C 

C«-«   GENERATE  PL0T3D  POINTER  FILE 
C 

DO  20  N=1 .NPLOTS+1 

XLOC  -  N»2.»PI/NPL0TS 
SINY(1,N)  -  (0.0)»YSCALE  +  YPOSIT 
SINX(1.N)  -  (XLOC)-XSCALE  +  XPOSIT 
SINY(2.N)  -  (SIN(XLOC)). YSCALE  +  YPOSIT 
SINX(2.N)  -  (XLOC)«XSCALE  +  XPOSIT 
20  CONTINUE 

PLOC  =  PNLOC • 2. »P I /NPLOTS 
SINY(3.1)  -  (0.0)»YSCALE  +  YPOSIT 
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SINX(3,1)  -  (PLOC)»XSCALE  +  XPOSIT 
SINY(3,2)  -  (SIN(PLOC))»YSCALE  +  YPOSIT 
SINX(3,2)  -  (PLOC)«XSCALE  +  XPOSIT 
C 

WRITE(50)  GDIM, NPLOTS 

KRITE(50)((SINX(I.J).  1-1 ,3). J-1 .NPLOTS). 
+  ((SINY(I.J).  1-1 ,3). J-1, NPLOTS) 

C 

PNLOC  -  PNLOC  +  1 

I F ( PN LOC . EO  NPLOTS ) THEN 

PNLOC  -  PNLOC  +  1  -  NPLOTS 
END  IF 
10     CONTINUE 
C 

C*«»   GENERATE  DUM^Y  0  FILE 
C 

WRITE(60)  GDIM, NPLOTS 
WRITE(60)  GDIM. GDIM. GDIM. GDIM 
WRITE(60)  ((01(1. J).  1=1 .3). J-1 .NPLOTS), 
+  ((02(1. J).  1-1 .3). J-1 .NPLOTS). 

+  ((03(1. J).  1=1 .3). J-1 .NPLOTS). 

+  ((04(1, J).  1-1. 3). J-1. NPLOTS) 

C 

1  FORMAT  (1X) 

2  FORMAT  (1X.2I8.4F8.4) 
END 

/EOF 
NPLOTS:  PNLOCI:  XSCALE:  YSCALE:  XPOSIT:  YPOSIT: 
294     222      .08      .1       -1.0    -.5 
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5.   CPPLOT.TXT 


JOB, JN-SKYHAWK.T-30.   ERIC  PAGENKOPF . X4269 . 
ACCOUNT  ,  AO ,  US-  .  UP*. . 

•  . 
CFT.ON-A.OFF-S. 

•  . 

CSF , ACCESS . CSFACCT- . CSFPSWD- . DN-GETP3D , PDN-GETP3D , I D-RFR I CC . 
ACCESS  ,  DN-SENDWKS  .  PDN-SENDWKS  .  I  D-STTRDM  .  OWN-RFTRDM  . 


CSF, ACCESS, CSFACCT-. CSFPSWD-.DN-CPX1 .  PDN-SNKR05P, 
CSF , ACCESS . CSFACCT- . CSFPSWO , DN=CFX1 .  PDN-SNKR05F , 
CSF . ACCESS . CSFACCT= . CSFPSWO . DN=Q1 .  PDN-SNKR05Q , 

REWIND.DN-CPX1 . 
DN=CFX1 . 
DN=01 . 

O-CPX . DATATYPE-XYZ . DSSK I P-49 , NUMDS-2 . 
O-CFX . DATATYPE-XYZ . DSSK I P-49 . NUMDS-2 . 
0-Q,       DATATYPE-Q.       DSSK I P-49 . NUMDS-2 . 


ID-FLYNAVY. 
ID-FLYNAVY. 
ID-FLYNAVY. 


REWIND 
REWIND 
GETP3D 
GETP3D 
CETP3D 


I-CPX1 
I-CFX1 

1=01  . 
REWIND. DN=CPV 
REWIND.DNNCFX 
REWIND. DN-Q. 
ASSIGN.DN-CPX.A 
ASSIGN.DN=CFX.A 
ASSIGN.DN=Q,       A 


■FT50. 
=FT60. 
■FT70. 


LDR . MAP-PART . SET-ZERO . 


/EOF 


PROGRAM  CPPLOT 


* 


THIS  PROGRAM  READS  STORED.  TRUE  VALUES  OF  XOC , 
AND  PRINTS  THEM  OUT  TO  TAPE  06. 


CP  AND  CF, 


VARIABLES: 
NUMDS: 
NPLOT : 


NUMBER  OF  DATA  SETS  IN  THE  READ  FILE 
PLOT  NO.  OF  FIRST  DATA  SET  (DSSKIP  +  1) 


TAPE  ASSIGNMENTS: 


TAPE  05 

TAPE  06 

TAPE  50 

TAPE  60 

TAPE  70 


INPUT  DATA 
OUTPUT  DATA 
TRUE  CP  VALUES 
TRUE  CF  VALUES 
0  FILES 


* 


DIMENSION  CPX(3.49),CFX(3.49).CPY(3,49).CFY(3.49) 

DIMENSION  01(161 ,41) .02(161 ,41). 03(1 61.41 ) ,Q4( 161 ,41) 

DIMENSION  KODE(4).LINE(90) 

DATA    K00E/1H    , 1H+, 1HI , 1H*/ 
C 

C»«»      READ    INPUT   DATA  It    INITIALIZE 
C 

READ    (5.1) 

READ  (5,2)  NUMDS. NPLOT 

DO  9  1=1 .90 

LINE(I)  -  K0DE(1) 
9     CONTINUE 
C 

C».»      READ  k   SAVE   TRUE  VALUES 
C 
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DO  10  L1-1 .NUMOS 


11 


ie 

c 

1 

2 
3 

4 

c 


READ(50) 
READ(50) 

READ (60) 
READ (60) 

READ(70) 
READ(70) 
READ(70) 


+ 
+ 
+ 


CP0 
CP0 

K0  i 
DO 


-  (1. 

=  CP0 
■*  30.  < 

11  L-1 , 


IMAX.KMAX 
((CPX(I,J 
((CPY(I ,J 
IMAX.KMAX 
((CFXM.J). 
((CFY(l.j). 
IMAXO.KMAXQ 
AMINF.ALFAD 

01(1. J). 
v i 02  C I . J ) . 
((03(1. J) 
((04(1. J) 


1-1 , IMAX) ,J-1 .KMAX), 
1-1 . IMAX) ,J-1 .KMAX) 


1-1 
1-1 


IMAX) ,J> 
IMAX). J- 


1 .KMAX) 
1 .KMAX) 


.REYREF.TIME 
1-1 .IMAXQ) .J 
1-1 . IMAXQ) 
1-1 . IMAXQ) 
1-1 . IMAXQ) 


+  .2 

/  (  • 

CP0 

KMAX 


AMINF..2) 
•  AMINF.»2) 
4.5 


1 .KMAXQ) 
J-1 , KMAXQ) 
J-1 .KMAXQ) 
J-1 .KMAXQ) 
5  -  1  . 


K  -  30.  •  (CP0  -  CPY(3.L))  +  4 

K1  -  30.  •  (CP0  -CPY(2.L))  +  4. 

IF(K.LT.1)  K  -  1 

IF(K1 .LT. 1)  K1  -  1 

IF(K.GT.90)  K  -  90 

IF(K1  .GT.  90)  K1  -  90 

LINE(K0)  -  K0DE(3) 

LINE(K)  =  K0DE(2) 

LINE(KI)  -  K0DE(4) 

IF(L.EQ.1)THEN 
WRITE(6.»)'  PLOT    AOA      TIME 
WRITE(6,3)  NPLOT.ALFAD. TIME. AMI NF.REYREF 
WRITE(6.»)'  XOC     CFL     CFU     CPL 
WRITE(6.4)  CPX(1 

ELSE 
WRITE(6.4)  CPX(1 ,L).CFY(2 

END  IF 


MACH    REY  NO. 


CPU 


L).CFY(2.L).CFY(3.L).CPY(2.L).CPY(3.L).LINE 
L).CFY(3.L).CPY(2.L).CPY(3.L).LINE 


LINE(K1)-K0DE(1 
KODE(V 


LINE(K) 
CONTINUE 
NPLOT  -  NPLOT 
WRITE(6.1) 
CONTINUE 


+  1 


FORMAT  (1X) 

FORMAT  (1X.2I7) 

FORMAT  (1X. 14.3F9.4.F10.0) 

FORMAT  (1X.F6.4.4F8.4.90A1) 


END 
/EOF 
NUMDS:  NPLOT: 
2      50 
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6.   GETX.F 


PROGRAM  GETX 


THIS  PROGRAM  READS  COMBINED  (MULTIPLE  DATA  SET)  FILES 
OF  XYZ  AND  0  PLOTTING  DATA  AND  SEPARATES  THEM  INTO 
INDIVIDUAL  FILES  FOR  PLOT3D  SUBMITTAL. 


VARIABLES: 

CFXNAME 
CFONAME 
FXNAME: 
FQNAME : 
DSSKIP: 

DSREAD: 
INTVL: 


NAME  OF  COMBINED  XYZ  FILE 

NAME  OF  COMBINED  Q  FILE 
NAME  OF  INDIVIDUAL  X  FILES.  CHAR  VARIABLE 
NAME  OF  INDIVIDUAL  Q  FILES.  CHAR  VARIABLE 
NUMBER  OF  DATA  SETS  TO  SKIP  WHEN 
READING  COMBINED  FILE 
NUMBER  OF  DATA  SETS  TO  READ 
INTERVAL  BETWEEN  DATA  SETS  SENT  TO  FILE 


INTEGER  DSSKIP.DSREAD. COUNT 

CHARACTER  FXNAME* 1 2 . FQNAME* 1 2 .CFXNAME* 10 , CFONAME*  1 0 

X(250. 100) .Z(250. 100) 

01(250. 100). 02(250. 100). 03(250. 100) .04(250. 


C 

•  • 

C 


DIMENSION 
DIMENSION 


100) 


•  USER  INPUT: 


C 

c 


DATA  XSCALE.YSCALE/1 . .1 ./ 

DATA  XP0SIT.YP0SIT/-1 . .-1  ./ 

DATA  DSSKIP/10/.DSREAD/3/. INTVL/1/ 

CFXNAME  -  '38x. iris' 

CFONAME  -  '38q. iris' 


0PEN(UNIT-91  .  FILE^FXNAME.STATUS^OLD'  ,  FORM- '  B I  NARY '  ) 
OPENtUNIT-^.FILE^FQNAME.STATUS-'OLD'  .  FORM=  '  B I  NARY  '  ) 


FXNAME 
FONAME 
COUNT 


•x000. 
q000. 
0 


iris 
iris1 


10 
C 


IF  (DSSKIP. GT.0)  THEN 

DO  10  ISKIP  -  1 .DSSKIP 


READ(91) 
READ(91) 

READ(92) 
READ(92) 
READ(92) 


+ 
+ 
+ 


IMAX.KMAX 

((X(I.K).I 

((Z(I.K).I 

IMAX.KMAX 

A.B.C.D 

((QUI.  K) 

((02(1. K) 
((03(1. K) 


IMAX),Ka 
IMAX) .K> 


CONTINUE 


END  IF 


1 ,KMAX) 
1 .KMAX) 


1  =  1 
1-1 

1-1 


((04(1. K). 1=1 


IMAX) 
IMAX) 
IMAX) 
IMAX) 


,K=1 ,KMAX) 
,K=1 .KMAX) 
pK=1 ,KMAX) 
,K=1 .KMAX) 


DO  20  IREAD  -  1 .DSREAD 
COUNT  =  COUNT  +  1 
WRITE  (FXNAME(2:4).100)  IREAD 
WRITE  (FQNAME(2:4),100)  IREAD 
READ(91)  IMAX.KMAX 
READ(91)  ((X(  I  .K).I-1  .  IMAX)  .K-=1 
I-  ((Z(1.K).I  =  1.  IMAX)  ,K=1 


.KMAX) 
,KMAX) 
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IF  (COUNT. EQ. INTVL)  THEN 

OPEN ( UN  I T-1 , FI LE-FXNAME . FORM- ' B I  NARY ' ) 
DO  11  1-1 , IMAX 

DO  12  K-1 .KMAX 

X(I ,K)-X(I .K)«XSCALE+XPOSIT 
Z(I.K)-Z(I.K)*YSCALE+YPOSIT 
12  CONTINUE 

11  CONTINUE 

WRITE(1)  IMAX, KMAX 

WRITE(I)  ((X(I .K) . 1-1 . IMAX). K-1 .KMAX) . 
+  ((Z(I.K). 1-1 .IMAX). K-1 .KMAX) 

CLOSE(1) 
END  IF 

READ(92)  IMAX. KMAX 
READ(92)  A.B.C.D 

READ(92)  ((01 (I ,K). 1=1 . IMAX).K=1 .KMAX), 
+  ((02(1 ,K). 1-1 , IMAX) .K-1 .KMAX). 

+  ((03(1 ,K). 1-1 . IMAX). K-1 .KMAX). 

+  ((04(1 .K). 1-1 . IMAX). K-1 .KMAX) 

IF  (COUNT. EO. INTVL)  THEN 

OPEN ( UN  I T-2 . F I LE=FONAME , FORM- ' B I  NARY ' ) 
WRITE(2)  IMAX. KMAX 
WRITE(2)  A.B.C.D 

WRITE(2)  ((01(1 ,K), 1-1 .IMAX) ,K-1 .KMAX), 
+  ((02(1 ,K) . 1-1 . IMAX) , K-1 .KMAX) . 

+  ((03(1, K).I«1 .IMAX). K-1 .KMAX). 

+  ((04(1. K), 1-1 , IMAX). K-1 .KMAX) 

CL0SE(2) 
COUNT  -  0 
END  IF 
20    CONTINUE 
C 

CL0SE(91) 
CL0SE(92) 
100    FORMAT (13.3) 
STOP 
END 
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7.   GETCP.F 


PROGRAM  GETCP 


•  THIS  PROGRAM  ACCESSES  COMBINED  CP   OR  CF  FILES  ON  THE  • 

•  IRIS  ACCOUNT.  SEPARATES  THEM  INTO  INDIVIDUAL  PLOTTING  • 

•  FILES  AND  SCALES  THEM  TO  MEET  SPECIFIED  REQUIREMENTS.  • 

•  A  DUMMY  0  FILr  !S  ALSO  GENERATED  (Q.IRIS).  • 

•  • 

•  VARIABLES:  ♦ 

•  » 

•  CPNAME:  COMBINED  FILE  NAME  • 

•  FNAME:  OUTPUT  FILE  NAME.  CHARACTER  VARIABLE  • 
.            DSSKIP:  NUMBER  OF  DATASETS  IN  THE  COMBINED  FILE  • 

•  TO  SKIP  • 

•  DSFILE:  NUMBER  OF  DATASETS  IN  THE  COMBINED  FILE  • 

•  TO  SEPARATE  AND  FILE  • 

AXIS  SCALING  FACTOR  • 

CP/CF  SCALING  FACTOR  (RELATIVE  MAGNITUDE)  • 

X  POSITION  FACTOR  (SCREEN  LOCATION)  • 

Y  POSITION  FACTOR  (SCREEN  LOCATION)  . 


•  XSCALE 

•  YSCALE 

•  XPOS I T 

•  YPOSIT 


INTEGER  DSSKIP. DSFILE 
CHARACTER  FNAME* 1 1 . QNAME»6 , CFNAME* 1 6 
DIMENSION  X(3.75).Z(3,75).XS(3.75).ZS(3.75) 
DIMENSION  01 (3. 75). Q2( 3. 75). 03(3.75) .04(3. 75) 
DATA  Q1 .Q2.Q3. 04/906* 1 ./ 


USER  INPUTS: 


DATA  DSSKIP/250/.DSFILE/44/ 
DATA  XSCALE/. 5/. YSCALE/-. 65/ 
DATA  XPOS I T/-1 .6/. YPOSIT/. 5/ 
CFNAME  -  " snkr65p" 
FNAME  ■  'cp606. iris' 


C 
C 

C 

•  •  • 
C 


OPENtUNIT-ge.FILE-CFNAME.STATUS-'OLD"  .  FORM- '  B I  NARY ' ) 
QNAME  -  ' q .  i  r  i  3  ' 
SKIP  DSSKIP  FILES 


IF  (DSSKIP. GT  6)  THEN 

DO  16  1SKIP  -  1 .DSSKIP 

READ(90)  IMAX.KMAX 

READ(96)  ( (X( I, K). 1-1. IMAX), K-1, KMAX) 

+  ((Z(I .K). 1-1. IMAX) , K-1 .KMAX) 
16           CONTINUE 
END  IF 


C 

•  •  • 

C 


SEPARATE  AND  SCALE  DSFILE  FILES 

DO  46  IFILE  -  1 .DSFILE 

WRITE(FNAME(3:5).106)  IFILE 
0PEN(UNIT«1 .FILE-FNAME. FORM- 'BINARY ' ) 
READ (96)  IMAX.KMAX 

READ(96)  ((X( I .K) . 1-1 . IMAX) .K-1 .KMAX) . 
+  ((Z(I .K) , 1-1 , IMAX) , K-1 .KMAX) 
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DO  30  I  ■  1 , IMAX 

DO  20  K  -  1 ,KMAX 

XS(I ,K)-X(I ,K)«XSCALE  +  XPOSIT 
ZS(I,K)=Z(I,K)»YSCALE  +  YPOSIT 
20  CONTINUE 

30  CONTINUE 

WRITE(I)  IMAX.KMAX 

WRITE(1)  ((XS(I ,K) , 1-1 . IMAX) ,K=1 , KMAX), 
+  ((ZS(I ,K). 1-1 , IMAX) ,K-1 ,KMAX) 

CLOSE(1) 
40      CONTINUE 


C 


GENERATE  DUMMY  Q  FILE 


OPEN ( UN  I T-1 . F I LE-ONAME . FORM-  ' B I  NARY ' ) 
WRITE(1)  IMAX.KMAX 
WRITE(1)  IMAX. IMAX, IMAX, IMAX 
WRITE(1)  ((01  U  .*)  .  1-1  .  IMAX)  ,K=1  .KMAX), 
+  ((02(1 ,K) ,1-1 , IMAX) ,K=1 .KMAX) . 

+  ((03(1 .*) .1=1 . IMAX) ,K«1 .KMAX) , 

+  ((04(1 ,K), 1-1 , IMAX) ,K=1 , KMAX) 

CL0SE(1) 
C 

100     FORMAT(I3.3) 
C 

STOP 
END 
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8.   GETSIN.F 


PROGRAM  GETS  IN 


THIS  PROGRAM  READS  A  COMBINED  FILE  OF  AOA  POINTER 
DATA  (GENERATED  WITH  PLOT3D  FUNCTION  0)  AND  SEPARATES 
IT  INTO  INDIVIDUAL  PLOTTING  FILES. 


VARIABLES: 
CFNAME 
FNAME: 
DSSKIP 


NAME  OF  COMBINED  FILE 
NAME  OF  INDIVIDUAL  FILES,  CHAR  VARIABLE 
NUMBER  OF  DATA  SETS  TO  SKIP  WHEN 
READING  COMBINED  FILE 
DSFILE:  NUMBER  OF  DATA  SETS  TO  SEND  TO  INDIVIDUAL 
FILES 


C 

•  •  • 


C 

c 
c 
c 


INTEGER  DSSKIP. DSFILE 
CHARACTER  FNAME» 1 2 .CFNAME* 1 0 
DIMENSION  X(3.294).Z(3.294) 


USER  INPUT: 


DATA  DSSKIP/250/.DSFILE/44/ 
CFNAME  =  'snkr05s' 


OPEN(UNIT=90.FILE=CFNAME.STATUS='OLD' , FORM- ' B I  NARY  * ) 
FNAME  -  'si  n000 .iris' 


IF  (DSSKIP. GT.0)  THEN 


DO  10  I  SKIP  - 

READ(90) 

READ(90) 

+ 

10 

CONTINUE 

END  IF 

c 

1 .DSSKIP 

IMAX.KMAX 

((X( I ,J) . 1=1 . IMAX).J=1 .KMAX) 

((Z(I . J) . 1—1 . IMAX).J=1 .KWAX) 


DO  20  IFILE  =  1 .DSFILE 

WRITE  (FNAME(4:6) ,100)  IFILE 
OPEN(UNIT-1 . FILE=FNAME,FORM= 'BINARY') 
READ(90)  IMAX.KMAX 
READ(90)  ((X(I.J),I«1 
h  ((Z(I.J).I-I 

IMAX.KMAX 
((X(I.J).I-1 . IMAX).J-1 .KMAX). 

( (Z( I.J). 1=1 . IMAX) ,  J=1 ,KUAX) 


IMAX).J=1 .KMAX). 
IMAX) . J=1 .KMAX) 


WRITE(1) 
WRITE(1) 

20 

CLOSE(1) 
CONTINUE 

C 
100 

CLOSE(90) 
FORMAT (I J. 3) 

STOP 

END 
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C.     PROCEDURAL  DOCUMENTATION 

After  proper  initialization  (user-inputs)  as  described  above,  all  codes  or  JCL's  which 
produce  or  transfer  Plot3D  files  aie  submitted  to  the  Cray  utilizing  standard  CSUB 
commands.  Once  the  data  is  resident  on  the  IRIS  and  the  user  is  logged  on  and  in  the 
proper  directory  level,  the  following  procedures  may  be  utilized  for  graphics  generation. 

1.  Separate  the  combined  file  into  discrete  data  sets. 

Edit  "user  inputs"  section  of  GE  IX. F  (or  other  files  as  appropriate)  to  define  scal- 
ing (ultimate  location  on  the  display),  and  identify  the  combined  file  to  access  and 
those  data  sets  to  be  separated: 

vi  getx.f 

Compile  the  editted  program: 

f77  getx.f 
Run  the  compiled  version  by  entering  the  filename  (a. out  by  default). 

a. out 

Discrete  data  sets  will  appear  on  the  account  with  names  corresponding  to  those 
in  the  Plot3D  initialization  files.    Check,  by  using  the  Is  command. 

2.  Use  Plot3D  to  generate  .gra  files. 

Use  the  mv  commands  to  give  the  proper  initialization  file  the  filename 
"plot3dini.com".  (The  IRIS  will  not  save  multiple  versions  of  files  with  the  same 
filename.  Instead,  older  versions  of  the  file  will  be  deleted.  In  order  to  avoid  in- 
advertent deletion  of  initialization  files,  multiple  mv  commands  may  be  required.) 

mv  xini.com  plot3dini.com 

If  processing  flow  field  files  (X  and  Q),  edit  the  initialization  fie  to  identify  desired 
function  and  number  of  contours.  This  involves  changing  only  two  lines  at  the  top 
of  the  file. 

vi  plot3dini.com 

Create  .gra  files. 

plot3x 
Rename  the  resultant  combined  graphics  file. 

mv  qOOl.gra  somefilnm.gra 

After  creation  of  as  many  .gra  files  as  required  from  the  separated  files,  clean  up 
the  account  using  the  following  wildcards. 

rm  q-.iris 

rm  x*.iris 
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3.    Interactive  viewing. 

Run  the  Graphics  Animation  System  software. 

gas 
Input  .gra  files  to  GAS  by  selecting  on  the  main  menu: 

ARCGRAPH  tile  input 
On  the  submenu,  select: 

load  entire  tile 

In  response  to  prompts,  enter  the  sequential  object  number  and  file  name 
("someilnm.gra").  For  color  plots,  press  RETURN  when  prompted  lor  the 
colormap. 

View  the  data  by  selecting  on  the  mam  menu: 

view  data 

Utilize  the  mouse  to  manipulate  the  display. 
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APPENDIX  B.     SANKAR  NAVIE      >TOKES  SOLVER 


(JCL) 

/EOF 


PROGRAM  MAIN 


SANKAR  NAVIER  STOKES  CODE:  SOLVES  TWO  DIMENSIONAL  FLOWS  PAST 
ARBITRARY  GEOMETRIES  USING  ADI  PROCEDURE.  THIS  VERSION  SAVES 
MULTIPLE  PLOT3D  DATA  SETS  IN  SINGLE  FILES. 


VARIABLES: 
TITLE 
I  MAX 

KMAX 
ITEL 
ITEU 
DT: 


(6C  CHARACTERS  MAX.) 
NUMBER  OF  X  COORD  LOCATIONS 
NUMBER  OF  Y  COORD  LOCATIONS 
GRID  POINT  AT  LOWER  TRAILING  EDGE 
GRID  POINT  AT  UPPER  TRAILING  EDGE 
SIZE  OF  TIME  STEP 


WW:  EXPLICIT  ARTIFICIAL  VISCOSITY  TERM 

ALFA:  MEAN  ANGLE  OF  ATTACK  (DEGREES) 

ALFA1 :  AMPLITUDE  OF  OSCILLATION  (DEGREES) 

ALFAI:  ANGLE  BELOW  WHICH  MODIFIED  TURBULENCE  MODEL  USED  TO 

COMPUTE  EDDY  VISCOSITY 
REDFRE:  REDUCED  FREOUENCY 
AM INF:  FREE  STREAM  MACH  NUMBER 

REYREF:  REYNOLDS  NUMBER  (MILLIONS;  NEG.  -  INVISCID  FLOW) 
DNMIN:  DISTANCE  OF  FIRST  POINT  OFF  OF  WALL 
TSTART:  TIME  FOR  CALCULATIONS  TO  START  ON  PRESENT  RUN 
FORMAT:  OLDSLN  FORMAT ; 3 . 0-PLOT3D  FILES, -3. 6-0  MATRICES 
RSTRT:  TRUE  =  STORED  DATA  READ  TO  CONTINUE  ITERATION 
PITCH:  TRUE  -  AIRFOL  AOA  OSCILLATES 
PLUNGE:  TRUE  =  AIRFOIL  OSCILLATES  VERTICALLY 
FNU:  NUMBER  OF  UPPER  SURFACE  DEFINITION  POINTS 
FNL:  NUMBER  OF  LOWER  SURFACE  DEFINITION  POINTS 
FSYM:  SYMMETRY  FLAG  (1  =  SYMMETRIC) 
X:  AIRFOIL  GEOMETRY  DEFINITION  POINTS 
Y:  AIRFOIL  GEOMETRY  DEFINITION  POINTS 


CSTP 
CPLT 
NSTP 
PSTP 


NUMBER  OF  STEPS  COMPLETED 

NUMBER  OF  DYNAMIC  PLOTS  COMPLETED 

NUMBER  OF  TIME  STEPS  TO  BE  COMPLETED  ON  PRESENT  RUN 

NUMBER  OF  TIME  STEPS  BETWEEN  PLOT  DATA  OUTPUT 


TAPE  DEFINIT 

TAPE  05 

TAPE  06 

TAPE  07 

TAPE  20 

TAPE  21 

TAPE  50 

TAPE  60 

TAPE  90 

TAPE  91 

TAPE  92 

TAPE  93 


IONS: 

INPUT  DATA 

OUTPUT  DATA 

FLOW  FIELD  INPUT  DATA  FOR  RESTARTS 

PL0T3D  XYZ  FILE  STORAGE 

PL0T3D  0  FILE  STORAGE 

PL0T3D  CP  XYZ  FILE  STORAGE 

PL0T3D  CF  XYZ  FILE  STORAGE 

PREVIOUS  RUN  X  FILE  STORAGE 

PREVIOUS  RUN  Q  FILE  STORAGE 

PREVIOUS  RUN  CP  FILE  STORAGE 

PREVIOUS  RUN  CF  FILE  STORAGE 


INTEGER  PSTP. CSTP, CPLT 
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c 

C*  »  • 

c 


c 

c*». 

c 


c 
c 


c 

c 
c 


c 

c..» 

c 


COMMON/SURF/PSUR  (161) 

COMMON/FIX/OMEGA 

COMMON/MUTUR/CMU (161,41) 

C0MM0N/SKINCF/CF(161) 

COMMON/GR I D 1 /X ( 1  6 1  ,41).Z(161,41) 

COMMON/PAR/GAMMA. REYREF. ALFA, ALFA1 . REDFRE . AMI NF , ALFAI 

COMMON/DGRID/DT. IMAX.KMAX, I  TEL. ITEU 

COMMON/GR  I  D/YACOB  (161.41) 

COMMON/DAMP/WW . WW  I 

COMMON/F LOW/01 (161  .41 ) .02(161  .41 ) .Q3(  161  .41 ) .04(161 .41) 

DIMENSION  TITLE(15) ,TYTLE(15) . ALPHA(96) ,CPTH(97 . 96) .XTH(97) 

DIMENSION  XR( 161 .49) ,ZR(161 .49) .Q1R( 1 61 .41) ,Q2R(161 .41) 

DIMENSION  Q3R(161 .41) ,04R(161 .41 ) 

COMMON/ LOG  I C/RSTRT . P I TCH . PLUNGE 

LOGICAL  RSTRT. PITCH, PLUNGE 

READ  INPUT  DATA 

READ  (5.5) 

READ  (5.1)  TITLE 

READ  (5.5) 

READ  (5.2)  IMAX.KMAX.DT.WW.ALFA.ALFA1 .ALFAI . REDFRE . AMINF 

READ  (5.5) 

READ  (5.3)  ITEL. ITEU. REYREF.DNMIN .TSTART . FORMAT. RSTRT. PITCH. PLUNGE 

READ  (5.5) 

READ  (5,4)  CSTP.CPLT.NSTP.PSTP 

INITIALIZE 

GAMMA  =1.4 

PI  =  ATAN(1 . )»4. 

SET  ALFAD  FOR  STEADY  STATE  PL0T3D  OUTPUT 

ALFAD  =  ALFA 

FORCE  DT  TO  BE  EOUAL  TO  UNITY  FOR  STEADY  FLOW  PROBLEMS 

THIS  INVOKES  THE  RELAXATION  MODE 

IF(REDFRE. LE. 0.001)  DT  =  1.0 

REYREF  =  REYREF  •  1 . E+06 

NITN  -  CSTP  +  NSTP 

NPLOTS  =  CPLT 

CPLT  «  CPLT  +  1 

CSTP  =  CSTP  +  1 

DENSITY  NORMALISED  WITH  RESPECT  TO  ROINF 

VELOCITIES  NORMALISED  WITH  RESPECT  TO  AINF 

TOTAL  ENERGY  NORMALISED  WITH  RESPECT  TO  (ROINF.AINF.AINF) 

TOTEN=AMINF. AMINF. 0.5+1 ./( GAMMA •( GAMMA- 1 . )) 

ALFA  =  ALFA  •  ATAN(1.)  /  45. 

ALMEAN  =  ALFA 

ALFAI  =  ALFAI  •  ATAN(1 

ALFA1  =  ALFA1  •  ATAN(1 

ALFAS  =  ALMEAN  -  ALFA1 

WRITE(6)'  PLOT   ITN 

UINF  =  AMINF 

VINF  =0.0 

00  7  1  =  1  . I  MAX 

DO  7  K«=1  .KMAX 

01 ( I ,K)=1 . 

02(1 .K)=UINF 

03(1 .K)=VINF 

Q4(I .K)=TOTEN 
CONTINUE 

INPUT  STORED  FLOW  FIELD  DATA 

IF(RSTRT)THEN 

I F( FORMAT. LT.0.0)THEN 

READ(7)  TIME. 01 .02.03.04 
ELSE 

READ(7)  IMAX.KMAX 


) 


45. 
45. 


TIME 


AOA 


CL 


CD 


CM' 
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READ(7) 
READ(7) 


AMI NF.ALFAD.REYREF, TIME 
((01(1. J),  1-1  .  IMAX).J=1  .KMAX) 


+ 
+ 
+ 


END  IF 
IF(NPL0TS 
DO  9 


((02(1. J). 
((03(1. J). 
((04(1. J). 


IMAX)  .  J- 
IMAX),J« 
IMAX)  .  J« 


1 .KMAX) , 
1 .KMAX) , 
1 ,KMAX) 


GT.0)THEN 
K1=1 .NPLOTS 


READ(90) 
READ(90) 

WRITE(20) 
WRITE(20) 


CONTINUE 
DO  10  K2=1 


IMAXR.KMAXR 
((XR(I.J).  I. 
((ZR(I.J).  I> 
IMAXR.KMAXR 
((XR(I.J).  I. 
((ZR(I.J).  I" 


, IMAXR) 
,  IMAXR) 

IMAXR) 
IMAXR) 


=  1  .KMAXR)  , 
«1  , KMAXR) 


J=1 , KMAXR) 
J-1 .KMAXR) 


NPLOTS 


+ 
+ 
+ 


+ 
+ 


READ(91) 
READ(91) 
READ(91) 


WRITE(21) 
WRITE(21) 
WRITE(21) 


IMAXR.KMAXR 

AM  I NFR . ALFADR , REYREFR . 
((01R(I.J).  1=1. IMAXR) 
((Q2R(I,J).  1-1, IMAXR) 
((Q3R(I.J).  1=1. IMAXR) 
((Q4R(I,J),  1=1. IMAXR) 
IMAXR.KMAXR 
AM I NFR, ALFADR 


TIMER 

. J=1 .KMAXR) , 
.J=1 .KMAXR). 
,J=1 .KMAXR), 
,J=1 .KMAXR) 


10 


CONTINUE 
DO  11  K3=1 


((Q1R(I.J) 

((02R(I.J) 
((Q3R(I.J) 
((04R(I.J). 


REYREFR, TIMER 
1-1 .IMAXR) ,J-1 .KMAXR). 
1-1. IMAXR). J-1 .KMAXR), 
1-1. IMAXR) ,J«1 .KMAXR), 
1=1 .IMAXR) ,J-1 .KMAXR) 


NPLOTS 


11 


READ (92)   IMAXR.KMAXR 

READ(92)   ((XR(I.J).  1-1 . IMAXR) , J-1 .KMAXR) , 

((ZR(I.J).  1=1 .IMAXR). J-1 .KMAXR) 
WRITE(50)  IMAXR.KMAXR 
WRITE(50)  ((XR(I.J).  1-1. IMAXR). J=1 , KMAXR). 

((ZR(I.J),  1=1 .IMAXR). J-1 .KMAXR) 


CONTINUE 

DO  12  K4-1 .NPLOTS 

READ(93)   IMAXR.KMAXR 

READ(93) 


WRITE(60) 
WRITE(60) 


((XR(I.J),  1=1 

((ZR(I.J).  1-1 

IMAXR.KMAXR 

((XR(I.J).  1=1 

((ZR(I.J).  1-1 


IMAXR), J=1 
IMAXR) , J-1 


KMAXR) 
KMAXR) 


IMAXR) 
IMAXR) 


12 


C 

C«  •  • 
C 

c 


J=1 .KMAXR) 
J-1  .KMAXR) 


15 


C 
C»  •  • 

c 


CONTINUE 
END  IF 
END  IF 

IF(TSTART.GE.0. )  TIME  -  TSTART 
IF( .NOT.(RSTRT))  TIME  -  0. 

GENERATE  COMPUTATIONAL  GRID.  DEFINE  SURFACE  COORD  O  0  AOA. 
ROTATE  TO  INITIAL  AOA  AND  COMPUTE  METRICS 

CALL  AIRFOL( IMAX. KMAX. ITEL. ITEU) 
IF  (REYREF.GT.0)  CALL  CLUSTR(DNMIN) 
TEDGE  -  X(ITEL  +  48,1) 
DO  15  I  -  1 .97 

XTH(I)  -  X(I+ITEL-1 .1)-TEDGE 
CONTINUE 

CALL  ROTGRID(X,Z, IMAX.KMAX , ALFAS) 
CALL  METRIC 

ITERATIVE  LOOP 

DO  1000  ITN-1 ,NSTP 
TIME  =  TIME  +  DT 
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<>••   ROTATE  GRID  TO  NEW  ANGLE  OR  SET  TO  CORRECT  ANGLE  FOR  RESTARTS 
C 

IF  (PITCH)  THEN 

OMEGA  -  2.  .  REDFRE  »AMINF.SIN(REDFRE  •  2.  •  TIME  •  AMINF) 
1       .ALFA1 

ALOLD  -  ALMEAN  -  ALFA1  •  COS(2.  •  REDFRE  •  AMINF  • 
1       (TIME  -  DT)) 

ALFA  -  ALMEAN  -  ALFA1  .  COS(REDFRE  •  2.  •  TIME  •  AMINF) 
ALFAD  =  ALFA  •  45 .  /  ATAN(1.) 
DALFA  =  ALFA  -  ALOLD 

IF(RSTRT.AND. ITN.EQ. 1 )  DALFA  =  ALFA  -  2.»ALFAS 
CALL  R0TGR1D(X.Z. IMAX.KMAX. DALFA) 
CALL  METRIC 
END  IF 
IF  (PLUNGE)  THEN 

OMEGA  =  2.  •  REDFRE  •  AMINF 
ALFA  -  ALMEAN 
END  IF 
C 

C»«»   COMPUTE  THE  SOLUTION  BY  AD  I  TECHNIQUE 
C 

CALL  SLPS(ITN. OMEGA. DALFA) 
C 

C...   APPLY  BOUNDARY  CONDITIONS 
C 

CALL  WALLBC 
C 

C»».   GENERATE  PLOT3D  FILES 
C 

IF(CSTP.EQ. (CPLT.PSTP))  THEN 
WRITE(20)  IMAX,  KMAX 

WRITE(20)  ((X(I,J).  1=1. IMAX).  J=1.KMAX). 
1  ( ( 2 ( I . J ) .  1=1. IMAX).  J-1.KMAX) 

WRITE(21)  IMAX,  KMAX 

WRITE(21)  AMINF.  ALFAD.  REYREF.  TIME 
WRITE(21)  ((Q1(I.J).  1-1. IMAX).  J-1.KMAX). 

1  ((Q2(I.J).  1=1. IMAX).  J-1.KMAX), 

2  ((Q3(I. J).  1-1. IMAX),  J-1.KMAX), 

3  ((04(1. J),  1=1. IMAX).  J- 1. KMAX) 
C 

C«««   GENERATE  PERFORMANCE  COEFFICIENTS 
C 

CALL  LOAD(PSUR.CL.CD.CM.ALFAS) 
AOA  =  ALFA* 180. /PI 

WR I T  E ( 6 , 6 )  CPLT . CSTP . T I  ME . AOA , C  L . CD . CM 
CALL  CPPL0T(ITEL.ITEU.AMINF,X(1 . 1 ) . 2( 1 .1) ,PSUR) 
CPLT  =  CPLT  +  1 
END  IF 
CSTP  =  CSTP  +  1 


1000 
c 

1 

CONTINUE 

FORMAT  (1X.15A4) 

2 

FORMAT  (1X.2I7.7F8.4) 

3 

FORMAT  (1X.2I7.4F8.5.3L7) 

4 

FORMAT  (1X.4I7) 

5 

FORMAT  (1X) 

6 

FORMAT  (1X, I3.I6.F8.3.F9.5.3F8 

C 

4) 


SUBROUTINE  AMAT1(K. IMX1 , XIX , XIZ ,XIT) 

C0MMON/FLOW/01(161 . 41 ) ,Q2( 1 61 . 41 ) ,Q3(  161  . 41 ) ,04( 1 61 . 41 ) 

C0MM0N/PERTR/DQ1(161 . 41 ) .DQ2( 161 .41 ) ,DQ3(  1  61 ,41).D04(161 .41) 

COMMON/AM/A (4, 4, 161) 

COMMON/PAR/GAMMA, REYREF. ALFA. ALFA 1 .REDFRE .AMINF, ALFAI 
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DIMENSION  XIT(161 
REAL  K0.K1 ,K2 


41),XIX(161 ,41),XIZ(161 ,41) 


C 

c 


AMAT1    COMPUTES   THE  COEFFICIENT  MATRIX   DE/DQ  DURING   XI    SWEEP 


GM1 

DO    1000 

K0 

K1 

K2 

U 

w 

EBYR 
PHI2 
THETA 
A(1.1.I) 
A(1.2.I) 
A(1,3.I) 
A(1.4,I) 
A(2.1.I) 
A(2.2.I) 
A(2.3.I) 
A(2,4.I) 
A(3.1 .1) 
A(3.2.I) 
A(3,3.I) 
A(3.4.I) 
A(4.1.I) 
A(4.2.I) 
A(4,3.I) 
A(4.4.I) 
1000  CONTINUE 
RETURN 
END 


GAMMA  -  1 . 

2  .  IMX1 

XIT(I.K) 

XIX(I.K) 

XIZ(I.K) 

Q2(I.K)  /  01(1. K) 

03(1. K)  /  01(1. K) 

04(1. K)  /  01(1. K) 


0.5  •  gmi  •  (u  •  u  +  w  •  w; 

) 

K1  •  U  +  K2  •  W 

K0 

K1 

K2 

0 

K1  •  PHI2  -  U  •  THETA 

K0  +  THETA  -  K1  •  (GM1  -  1 

) 

•  U 

K2  •  U  -  GM1  •  K1  •  W 

K1  •  GM1 

K2  •  PHI2  -  W  •  THETA 

K1  •  W  -  K2  •  GM1  •  U 

K0  +  THETA  -  K2  •  (GM1  -  1 . 

) 

•  W 

K2  •  GM1 

THETA  •  (2.  •  PHI2  -  GAMMA 

* 

EBYR) 

K1  «  (GAMMA  •  EBYR  -  PHI2) 

- 

GM1  •  U  • 

THETA 

K2  •  (GAMMA  •  EBYR  -  PHI2) 

- 

GM1  •  W  • 

THETA 

K0  +  GAMMA  •  THETA 

c 
c« 

c 


SUBROUTINE  AMAT2(I .KMX1 , ZETAX.ZETAZ.ZETAT) 

COMMON/ FI0W/Q1 (161 .41 ) ,02(161 .41 ) .03(161 ,41 ) ,04(161 .41 ) 

C0MM0N/PERTR/DQ1(161 ,41 ) .002(161 ,41 ) ,003(161 ,41). 004(161 .41) 

COMMON/ AM/ A( 4. 4. 161 ) 

COMMON/PAR/GAMMA . REYREF . ALFA . ALFA  1 . REDFRE . AMI NF , ALFAI 

DIMENSION  ZETAX(161 . 41 ) ,ZETAZ( 161 .41 ) ,ZETAT( 161 ,41) 

REAL  K0.K1 ,K2 


c 

c  •  •  • 

c 

AMAT2  COMPUTES  THE  COEFFICIENT  MATRIX  DF/D 

GM1 

= 

GAMMA  -  1 . 

DO  1000  K 

= 

2  .  KMX1 

K0 

= 

ZETAT(I.K) 

K1 

■ 

ZETAX(I.K) 

K2 

= 

ZETAZ(I.K) 

U 

= 

02(1. K)  /  01(1. K) 

w 

= 

03(1. K)  /  01(1, K) 

EBYR 

= 

04(1. K)  /  01(1, K) 

PHI2 

s 

0.5  •  GM1  •  (U  •  U  +  W  •  W) 

THETA 

■ 

K1  •  U  +  K2  •  W 

A(1.1.K) 

= 

K0 

A(1.2.K) 

- 

K1 

A(1.3,K) 

= 

K2 

A(1 ,4.K) 

■ 

0 

A(2.1,K) 

a 

K1  •  PHI2  -U  •  THETA 

A(2.2.K) 

» 

K0  +  THETA  -  K1  •  (GM1-1.)  •  U 

A(2.3.K) 

= 

K2  •  U  -  GM1  •  K1  •  W 

A(2.4.K) 

= 

K1  •  GM1 

A(3,1 ,K) 

= 

K2  •  PHI 2  -  W  •  THETA 

A(3,2,K) 

= 

K1  •  W  -  K2  •  GMI  •  U 

A(3.3.K) 

■ 

K0  +  THETA  -K2  •  (GM1-1 . )  •  W 
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A(3.4.K) 
A(4,1 .K) 
A(4.2.k) 
A(4.3.K) 
A(4,4,K) 
1600  CONTINUE 
RETURN 
END 


K2  •  GM1 

THETA  •  (2.  •  PHI2  - 

-  GAMMA 

•  EBYR) 

K1  •  (GAMMA  •  EBYR  - 

-  PHI2) 

-  GM1  •  U  •  THETA 

K2  •  (GAMMA  •  EBYR  - 

-  PHI2) 

-  GM1  •  W  •  THETA 

K0  +  GAMMA  •  THETA 

c 


SUBROUTINE  SLPS(ITN. OMEGA. DALFA) 

COMMON/F LOW/01 (  1  61  . 41 ) . Q2( 1 61 . 41 ) . Q3( 1 61  . 41 ) , Q4( 1 61  . 41 ) 

COMMON/PERTR/DQ1(161 ,41 ),DQ2(161 .41), 003(161 ,41) ,DQ4(161 .41) 

COMMON /AM/A (4, 4. 161 ) 

COMMON/TRID/DD(4,4.161 , 41 ) ,MM(4 , 4 . 161 . 41 ) . EE(4.4, 1 61 .41 ) 
1  ,GG(4.  161  ,41 ) 

COMMON/PAR/GAMMA . REYREF , ALFA . ALFA1 . REDFRE . AMINF , ALFA I 

COMMON/DGR I D/DT , I  MAX . KMAX , I  TEL . I TEU 

CCMMON/GR  I D/YACOB  (161.41) 

COMMON/DAMP/WW , WW  I 

COMMON/MTRIX/  XIX(161,41),XIZ(161 . 41 ) . ZETAX( 1 61 ,41), 
+ZETAZ(161 ,41) 
1  ,XIT(161 ,41) ,ZETAT(161 .41) 

REAL  MM 

DIMENSION  DELTAT(161 .41) 
C 

C-«-   SUBROUTINE  SLPS  DOES  THE  BULK  OF  THE  WORK   FOR  THE  AC:  ALGORITHM. 
C«««   IT  CALLS  FLUX  AND  COMPUTES  RIGHT  HAND  SIDE  DURING  THE  TWO  SWEEPS. 
C««»   ASSEMBLES  THE  COEFFICIENT  MARICES.  ADDS  IMPLICIT  AND  EXPLICIT 
C««»   DISSIPATION  AND  CALLS  THE  TRIDIAGONAL  SOLVER  TO  OBTAIN  THE  FINAL 
C».«   SOLUTION. 
CMM.'SET  VALUE  OF  IMPLICIT  DAMPING  COEFFICIENT 

WWI  =  20.0  •  WW 

IM1  -  IMAX  -  1 

IM2  -  IMAX  -  2 

KM1  =  KMAX  -  1 

KM2  »  KMAX  -  2 

IF(ITN.EQ.I)  THEN 
C! ! ! ! ! LOCAL  TIME  STEP  OPTION  FOR  STEADY  STATE  CALCULATIONS 

IF(REDFRE.LT. 0.001)  THEN 

DO  777  K  -  2  ,  KMAX  -  1 

DO  777  1=2.  IMAX  -  1 

DELTAT(I.K)  =  0.5  /  (1.  +  SORT(ABS(YACOB( I ,K)) ) ) 

777  CONTINUE 
ELSE 

DO  778  K  =  2  .  KMAX  -  1 
DO  778  1=2.  IMAX  -  1 

778  DELTAT(I  ,K)  =  1 .0 
END  IF 

END  IF 
C 

C*»«   THE  DISSIPATION  TERMS  ARE  CONSTRUCTED  AND  STORED  IN  THE  ARRAYS  DQ1 . 
C«««   DQ2.DQ3  AND  DQ4. 
C 

CALL  DISSIP 
C 
C»«»   THE  RIGHT  HAND  SIDE  AT  KNOWN  TIME  LEVEL  IS  NOW  COMPUTED  AND  ADDED 

CALL  RESI (OMEGA) 
C 
C»*»       IF  VISCOUS  FLOW  IS  COMPUTED  THE  VISCOUS  TERMS  ARE  ADDED  TO  DQ1  ETC.  HERE. 


C 

C< 

c 


IF(REYREF.GT0. )  CALL  STRESS( ITN. DALFA) 
I -SWEEP. 

DTH  =  DT  .  0.5 
DTW  =  DT  .  WWI 
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DO  3 
CALL 
DO  4 
DO  4 
DO  5 
EE(  1 1 
DD(I1 
CONTINUE 
CONTINUE 


K  -  2  .  KM1 

AMAT1(K. IMAX-1 .XIX.XIZ.XIT) 


II 
12 
I 


12. 1-1. K) 
12. 1-1 .K) 


1  .  4 

1  .  4 

2  .  IMAX  -  1 
A(I1 .12.1+1) 

-  A( 1 1 .12.1-1) 


DTH 
DTH 


DELTAT(I,K) 
DELTAT(I.K) 


C 

C»»»   IMPLICIT  DAMPING  ADDED  HERE. 

C 


990 
3 


1  .  4 

2  .  IMAX  -  1 

DTW  /  YACOB(I.K)  «  DELTAT(I.K) 
D0( 1 1 . 11 . 1-1 .K)  -  DT1  •  YACOB(I-I.K) 
EE(I1 .11 .I-1.K)  -  DT1  .  YACOB(I+1.K) 
1 .  +  2.  •  DTW  •  DELTAT(I.K) 

2  ,  IMAX  -  1 
DQ1(I ,K)  •  DELTAT(I ,K) 
D02(I.K)  •  DELTAT(I.K) 
DQ3(I ,K)  •  DELTAT(I ,K) 
D04(I.K)  •  DELTAT(I.K) 


C 

C 
C 


c 

C 

c 


15 


C 

C« 

C 


c 

C 

c 


DO  6  11 

DO  6  I 

DT1 

DD( 1 1  ,11.1-1 .K) 

EE(I1  .11.1-1 .K) 

MM(I1 .11.1-1 ,K) 

CONTINUE 

DO  990  I 

GG(1 . I-1.K) 

GG(2.I-1.K) 

GG(3. 1-1. K) 

GG(4.I-1.K) 

CONTINUE 

CONTINUE 

PERFORM  BLOCK  TRIDIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 

CALL  MATRX1  (IMAX.KMAX) 


DO  991  K 
DO  991  I 
D01(I ,K) 
DQ2(I.K) 
D03(I.K) 
D04(I.K) 
991  CONTINUE 


2  ,  KMAX  -  1 
2  ,  IM1 
GG(1 . 1-1 .K) 
GG(2. 1-1 .K) 
GG(3,I-1 ,K) 
GG(4,I-1 ,K) 


K-SWEEP  BEGINS  HERE. 

DO  13  I         =  2  ,  IM1 

CALL  AMAT2(I .KMAX-1 .2ETAX . ZETA2 . ZETAT) 


DO  15  II 

DO  15  12 

DO  15  K 

EE(I1  .I2.I.K-1) 

DD(I1  .12.1 ,K-1) 

CONTINUE 


1  .  4 

1  .  4 

2  .  KMAX  -  1 

A(I1 ,I2.K+1)»DTH  •  DELTAT(I.K) 
-A(I1 ,I2.K-1).DTH  •  DELTAT(I.K) 


SECOND  ORDER  DAMPING  ADDED  HERE. 


DO  16  11 

DO  16  K 

DT1 

DD(I1  .I1.I.K-1) 

EE(I1  .I1.I.K-1) 

16  MM(I1  .  H. I  .K-1) 
DO  17  K 

GG(1  .  I  .K-1) 
GG(2.I.K-l) 
GG(3. I .K-1) 
GG(4,I  ,K-1) 

17  CONTINUE 
13  CONTINUE 


1  .  4 

2  ,  KMAX  -  1 
DTW  /  YACOB(I.K) 
DD(I1 .11 .1 .K-1)  - 
EE(I1 .11.1 
1.  +  2.  • 
2  .  KMAX  - 
D01(I.K) 
D02(I.K) 
DQ3(I,K) 
D04(I.K) 


DELTAT(I.K) 
DT1  •  YACOB(I .K-1) 
K-1)  ■-  DT1  •  YACOB(I  ,K+1) 
DTW  •  DELTAT(l.K) 
1 


PERFORM  BLOCK  TRIDIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 
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CALL  MATRX2(IMAX.KMAX) 


DO  18  I 

-  2  .  IMAX  -  1 

DO  18  K 

-  2  ,  KM1 

D01(I .K) 

-  GG(1 .1 .K-1) 

DQ2(I.K) 

-  GG(2.I.K-1) 

D03(I .K) 

-  GG(3.I ,K-1) 

DQ4(I .K) 

-  GC(4.I .K-1) 

18 
C 
C»  •  • 

CONTINUE 

UPDATE  FLOW 

VARIABLES  AT  INTERIOR 

POINTS. 

967 

CONTINUE 

RMAX 

RUMAX 

RVMAX 

EMAX 

DO  995  K 

DO  19   I 

-  0. 
«=  0. 

-  0. 
=  0. 

=  2  .  KM1 

-  2  ,  IM1 

01(1. K) 

-  01(1 ,K)  +  D01(I 

,K)  •  YACOB(I 

•  K) 

Q2(I.K) 

-  02(1 .K)  +  DQ2(I 

,K)  •  YACOB(I 

•  K) 

03(1. K) 

=  03(1 ,K)  +  DQ3(I 

,K)  •  YACOB(I 

•  K) 

04(1 .K) 

-  04(1 .K)  +  DQ4(I 

,K)  •  YACOB(I 

•  K) 

19  CONTINUE 
CIMMDETE  MINE  WHERE  IN  FLOW  FIELD  DENSITY  IS  CHANGING  MOST  RAPIDLY 
DO  Jd5    I         =  2  .  IMAX  -  1 
IF  (RMAX.LT.ABS(DQI(I.K)-YACOB(I .K)))  THEN 
IR  -  I 

KR  -=  K 

END  IF 
RMAX 
RUMAX 
RVMAX 
EMAX 
CONTINUE 

IF((ITN-1)/100«100.EQ. (ITN-1))  WRITE  (6.3002) 

IF(ITN  .EO.  0)  WRITE  (6.3002) 
SELECT  INTERVAL  AT  WHICH  OUTPUT  OF  RESIDUALS  IS  DESIRED 

IF((ITN-1)/10»10.EQ.(ITN-1))  WRITE  (6.3001)  RMAX , RUMAX , RVMAX . 
1EMAX.IR.KR 
RETURN 

FORMAT (// . 4X . *  DRMAX ' .  1 1 X . ' DUMAX ' . 1 1 X . ' DVMAX ' . 1 1 X . ' DEMAX ' . 1 0X . 
1 • IR' .3X, "KR") 
3001  F0RMAT(4(E14.8.2X),2I5) 
END 


995 
C 
C 

C!  !  !  !  ! 
C 
C 

3002 


AMAX1 (RMAX.ABS(D01 (I ,K)  •  YACOB(I.K))) 
AMAX1 ( RUMAX. ABS(DQ2( I ,K)  •  YACOB(I.K))) 
AMAX1 ( RVMAX. ABS(DQ3( I .K)  •  YACOB(I.K))) 
AMAX1(EMAX.ABS(D04(I ,K)  .  YACOB(I.K))) 


C 
C< 

c 


SUBROUTINE  MATRX1 ( IMAX.KMAX) 

COMMON/TRID/DD(4.4.161  .  41  ) ,MM(4. 4 . 1 61 . 41 ) . EE(4. 4 , 1 61 . 41 ) , 
1GG(4.161 .41) 

C0MM0N/SCRAT/A(4.4.161) ,HH(4,4, 161 .41 ) ,C(4.5. 161 ) 

REAL  MM 

REAL  L11 . L21 . L31 . L41 . L22 . L32 . L42 . L33 . L43. L44 
2.L1I  ,L2I . L3I .L4I 


C 
C 
C 
C 


THIS  SUBROUTINE  PERFORMS  THE  BLOCK  TRIDIAGONAL  MATRIX  IVERSION  FOR 
AN  ENTIRE   PLANE  DURING  THE  XI-  SWEEP 

DO  1  11-1.  4 

DO  1  K  -  2  ,  KMAX  -  1 

AI  =  1 .  /  MM(1 .1 .1 ,K) 


GG(I1.1.K)  = 


HH(I1 .1 .1 
HH(I1 .2,1 
HH(I1 .3.1 
HH( 11 .4.1 
CONTINUE 


,K) 
,K) 

■  K) 

,K) 


GG ( 1 1 . 1 .K)  •  AI 

-  EE(I1.1 .1 ,K)  •  AI 

-  EE( 11.2.1  .K)  •  AI 

-  EE( 11.3.1  ,K)  •  AI 

-  EE(I1.4.1 ,K)  •  AI 
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DO  1000  I  - 
DO  5  II- 
DO  2  K  « 
C(II.I.K) 

1 

2 
3 
2  CONTINUE 


.  IMAX 
,  4 

.  KMAX 
GG ( 1 1 . I 


-  2 

-  1 
.K)  - 


DD(I1 .1 . I.K)  •  GG(1 .1-1 
DD(I1 .2. I.K)  •  GG(2.I-1 
DD(I1 ,3. I ,K)  •  GG(3. 1-1 
DD(I1 .4. I.K)  »  GG(4.I-1 


DO  5 
DO  5 
A(  I  1 


12  ■ 

K   ■ 
12. K). 


.  4 

.  KMAX 
(11.12. 


-  1 
I.K)  - 


1 

2 
3 


C(I1.I2+1.K)-  EE(I1 .12. I.K) 

CONTINUE 

DO  3  K  -  2 


DD(I1 .1 , I.K) 
DD(I1 .2. I.K) 
DD( 1 1 .3. I.K) 
DD(I1 .4, I.K) 


HH(1.I2, 
HH(2.I2, 
HH(3.I2, 
HH(4. 12. 


K) 
K) 
K) 
K) 


1-1. K) 
1-1 .K) 
1-1 .K) 
1-1 .K) 


L11 
L1I 
U12 
U13 
U14 
L21 
L31 
L41 
L22 
L2I 
U23 
U24 
L32 
L42 
L33 
L3I 
U34 
L43 
L44 
L4I 

CO 


,  KMAX  - 
,K) 


1 


L1I 
L1I 
L1I 


-  L21 


1 


A(1  .1 

1 .  /  L11 

A(1 ,2.K) 

A(1 .3.K) 

A(1 ,4,K) 

A(2.1 .K) 

A(3.1 .K) 

A(4.1 ,K) 

A(2.2.K) 

1 .  /  L22 

(A(2.3.K)  -  L21 

(A(2.4.K)  -  L21 

A(3.2.K)  -  L31 

A(4,2.K)  -  L41 

A(3.3.K)  -  L31 

1 .  /  L33 

(A(3.4.K)  -  L31 

A(4.3.K)  -  L41 

A(4,4.K)  -  L41 

1 .  /  L44 

K)  -  C(1 .1 ,K)  • 
(C(2.1 ,K) 
(C(3.1 ,K) 


»  U12 

•  U13 

•  U14 

•  U12 

•  U12 

•  U13 

•  U14 

•  U13 

•  U14 


*  L2I 

•  L2I 


-  L32 


U23 


L3I 


-  L32  •  U24) 

-  L42  •  U23 

-  L42  •  U24  -  L43  •  U34 


L1I 
C(2.1 ,K)  =  (C(2.1 ,K)  -  L21 
C(3.1.K)  »  (C(3.1 ,K)  -  L31 

-  L32 
C(4  1 ,K)  -  (C(4.1 ,K)  -  L41 


C  J.K) 
C  .K) 
C    .K) 

C(4.2.K) 

C(1.3,K) 
C(2.3.K) 
C(3.3.K) 

C(4.3,K) 

C(1  ,4.K) 
C(2.4,K) 
C(3.4.K) 

C(4.4.K) 

C(1 .5.K) 
C(2.5.K) 
C(3.5.K) 


C(1.2.K)  •  L1I 
(C(2.2.K)  -  L21 
(C(3,2,K)  -  L31 

-  L32 
(C(4.2.K)  -  L41 

C(1 .3.K)  •  L1 I 
(C(2.3,K)  -  L21 
(C(3.3,K)  -  L31 

-  L32 
(C(4.3.K)  -  L41 

C(1 .4.K)  •  L1I 

(C(2.4.K)  -  L21 
(C(3,4.K)  -  L31 

-  L32 
(C(4.4.K)  -  L41 


C(1.5.K) 
(C(2.5.K) 
(C(3.5.K) 


C(4.5.K)  =  (C(4.5.K)  - 


LI  I 

L21 
L31 
L32 
L41 


C(I.I.K)) 
C(1 .1 ,K) 
C(2.1 ,K)) 
C(1 .1 .K)  - 

-  L43  • 

C(1.2.K)) 
C(1.2.K) 
C(2.2.K)) 
C(1.2.K)  - 

-  L43  • 

C(1.3.K)) 

C(1 .3.K) 
C(2.3.K)) 
C(1.3.K)  - 

-  L43  • 

C(1.4.K)) 
C(1 .4.K) 
C(2.4.K)) 
C(1  ,4.K)  - 

-  L43  • 

C(1 .5.K)) 
C(1 .5.K) 
C(2.5.K)) 


L2I 

L3I 
L42  • 
C(3.1 

L2I 

L3I 
L42  • 
C(3.2 

L2I 


C(2( 
K)) 


C(2.2 
K))  • 


L3I 
L42  •  C(2.3 
C(3.3.K))  • 

L2I 

L3I 

L42  •  C(2,4, 
C(3.4.K))  ♦ 

L2I 

L3I 


C(1 ,5.K)  -  L42  •  C(2.5 
-  L43  •  C(3.5.K))  • 


■  K) 
L4I 


K) 

L4I 


K) 
L4I 


K) 

L4I 


K) 

L4I 
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Cf  3 . 1 .K) 

C(2.1 ,K) 


C(3.1  ,K) 
C(2.1.K) 


C(1 ,1 ,K)  -  C(1 .1 ,K)  - 


C(3.2.K) 


C(3.2.K) 
C(2.2,K) 


C(2.2,K 

C(1 ,2.K)  -  C(1 ,2.K)  - 


C(3.3.K) 
C(2.3.K) 

C(1 ,3.K) 

C(3.4.K) 
C(2.4.K) 

C(1 ,4.K) 

C(3.5.K) 
C(2.5.K) 

C(1 .5.K) 

3  CONTINUE 

DO  6  II 
DO  9  K 

9  GG(I1  .1 
DO  6  12 
DO  6  K 
HH(I1 .12, 

6  CONTINUE 
1000  CONTINUE 


C(3.3.K) 
C(2.3,K) 


C(1 ,3.K) 

C(3.4.K) 
C(2.4.K) 

C(1 ,4.K) 


C(3.5.K} 
C(2.5.K) 

C(1 ,5,K) 


U34 
U24 
U23 
U14 

ui: 

U3* 
U24 

U23 
U14 
U13 
U34 
U24 
U23 
U14 
U13 
U34 
U24 
U23 
U14 
U13 
U34 
U24 
U23 
U14 
U13 


C(4. 
C(4. 
C(3. 
C(4, 
C(3. 
C(4. 
C(4. 
C(3. 
C(4. 
C(3. 
C(4. 
C(4. 
C(3. 
C(4. 
C(3. 
C(4, 
C(4, 
C(3. 
C(4. 
C(3. 
C(4. 
C(4. 
C(3, 
C(4. 
C(3. 


1,K) 

1.K) 

I.K 

I.K) 

1.K) 

2.K) 

2.K) 

2.K) 

2.K) 

2.K) 

3.K) 

3.K) 

3.K) 

3.K) 

3.K) 

4.K) 

4.K) 

4.K) 

4.K) 

4.K) 

5.K) 

5.K) 

5.K) 

5.K) 

5.K) 


-  U12  •  C(2.1 .K) 


-  U12  •  C(2.2,K) 


-  U12  •  C(2.3.K) 


-  U12  •  C(2,4.K) 


-  U12  •  C(2.5.K) 


■1,4 

-  2  .  KMAX  -  1 
K)     ■  C(I1,1,K) 

■  1,4 

-  2  .  KMAX  -  1 
I.K)  -  C(I1,I2+1.K) 


C 
C 
C 


BACKWARD  SUBSTITUTION 


DO  7  I 

DO  7  II 

DO  7  K 

GC(I1 .I.K) 
1 
2 
3 
7  CONTINUE 

RETURN 

END 


IMAX  -  3.  1  .  -  1 

1  .  4 

2  .  KMAX  -  1 
GG(M  .I.K)  -  HH(I1 

-  HH(I1 
HH(I1 


I.K) 
I.K) 
I.K) 


-  HH(I1 .4. I.K)  • 


GG(1.I+1 ,K) 
GG(2.I+1 .K) 
GG(3.I+1 ,K) 
GG(4.I+1.K) 


C 

C< 
C 


c 
c 
c 
c 


SUBROUTINE  MATRX2( IMAX .KMAX) 

COMtf)N/TRID/DD(4.4.161 .41 ) .1*1(4.4. 161 .41 ) . EE(4.4. 161 .41). 
1GG(4.161  .41) 

COMMON/SCRAT/A(4.4.161).HH(4,4.161 . 41 ) ,C(4,5. 161  ) 

REAL  MM 

REAL  L11 .L21.L31.L41.L22.L32.L42.L33.L43.L44 
2.L1I.L2I.L3I.L4I 

THIS  SUBROUTINE  PERFORMS  THE  BLOCK  TRIDIAGONAL  MATRIX  IVERSION  FOR 
AN  ENTIRE  J-CONSTANT  PLANE  DURING  THE  2ETA-  SWEEP 

DO  1  II  -1,4 

DO  1  I  -  2  .  IMAX  -  1 

AI  «  1 .  /  MM(1. 1.1.1) 

GG(I1  .1.1)  -  GG ( 1 1 . 1 . 1 )  •  AI 

HH(I1.1  .1.1)  -  EE(I1. 1.1.1)  •  AI 
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HH(I1 .2. I .1)  -  EEC  II .2.1.1)  •  AI 

HH(I1 .3, I .1)  =  EE( 1 1 .3.1 .1)  •  AI 

HH(I1 .4.1 .1)  -  EE(I1 ,4. I . 1)  •  AI 

1  CONTINUE 

do  ieee  k  -  2  .  kmax  -  2 

DO  5     11-1,4 

DO  2     I  -  2  .  IMAX  -  1 

C(  1 1  .  1  .1)  -  GG(II.I.K)  -  OD(I1.1. 1, K)  •  GG(1.I.K-1) 

1  -  DD(I1  .2.I.K)  •  GG(2.I.K-1) 

2  -  DD(I1 .3.1 ,K)  •  GG(3.I,K-1) 

3  -  DD(11,4.I.K)  •  GG(4,I,K-1) 

2  CONTINUE 

DO  5    12  «=  1  .  4 

DO  5    I   -  2  ,  IMAX  -  1 

A( 1 1.12. 1)-  MI(I1.I2.I.K)  -  DD(I1  .1.I.K)  .  HH(1 . 12 . 1 .K-1 ) 

1  -  00(11  .2.1. K)  •  HH(2.I2.I.K-1) 

2  -  DDCH.3.I.K)  •  HH(3.I2,I.K-1) 

3  -  DD(I1.4.I.K)  .  HH(4.I2.I.K-1) 
C(I1 .12+1.1)-  EEC  1 1 . 12. I ,K) 

5  CONTINUE 

DO  3  I  -  2  .  IMAX  -  1 

L11  =  A( 1  . 1  . 1 ) 

L1I  -  1.  /  L11 

U12  -  A(1 .2.1)  •  L1 I 

U13  =  A(1 .3.1)  •  L1I 

U14  =  A(1 .4.1)  •  L1 I 

L21  -=  A(2.1.I) 

L31  «=  A(3.1  .1) 

L41  =  A(4.1 .1) 

L22  -  A(2.2.I)  -  L21  •  U12 

L2I  -  1  .  /  L22 

U23  -  (A(2.3.I)  -  L21  •  U13)  •  L2I 

U24  -  (A(2.4.I)  -  L21  •  U14)  «  L2I 

L32  -  A(3.2.I)  -  L31  •  U12 

L42  -  A(4.2. I)  -  L41  •  U12 

L33  -  A(3,3,I)  -  L31  •  U13  -  L32  •  U23 

L3I  -  1.  /  L33 

U34  -  (A(3.4.I)  -  L31  •  U14  -  L32  •  U24)  •  L3I 

L43  «=  A(4.3,I)  -  L41  •  U13  -  L42  •  U23 

L44  -  A(4.4.I)  -  L41  •  U14  -  L42  •  U24  -  L43  •  U34 

L4I  ■  1 .  /  L44 

C(1 ,1  .1)  =  C(1,1 . I)  •  L1I 

CC2.1.I)  =  (C(2.1.I)  -  L21  •  C(1.1.I))  •  L2I 

C(3.1.I)  =  (C(3.1.I)  -  L31  »  C(1.1.I) 

-  L32  *  C(2.1  .1))  •  L3I 
C(4.1,I)  *  (C(4.1.I)  -  L41  •  C(1.1.I)  -  L42  »  C(2.1.I) 

-  L43  •  C(3.1  .1))  .  L4I 
CO  .2.1)  «  C(1  .2.  I )  •  L1I 

C(2.2.I)  =  (C(2.2.I)  -  L21  •  C(1.2.I))  •  L2I 

C(3.2.I)  =  (C(3.2.I)  -  L31  •  C( 1 .2.1) 

-  L32  •  CC2.2.I))  •  L3I 
C(4.2.I)  =  (C(4.2,I)  -  L41  .  C(1 .2. I)  -  142  •  C(2.2.I) 

-  L43  •  C(3.2.I))  •  L4I 
C(1 .3.1)  =  C(1.3.I)  •  L1I 

C(2.3.1)  -  (C(2.3.I)  -  L21  •  C(1.3.I))  •  L2I 

C(3.3.I)  =  (C(3.3.I)  -  L31  •  C(1.3.I) 

-  L32  •  C(2.3.I))  •  L3I 
C(4.3,I)  -  (C(4.3.I)  -  L41  .  C(1.3.I)  -  L42  •  C(2.3.I) 

-  L43  •  C(3.3.I))  •  L4I 
C(1 .4.1)  -  C(1.4.I)  •  L1I 

C(2.4.I)  =  (C(2.4.I)  -  L21  »  C(1.4.I))  •  L2I 

C(3.4.I)  «  (C(3.4.I)  -  L31  •  C(1.4.I) 

-  L32  •  C(2.4.I))  •  L3I 
C(4.4.I)  -  (C(4,4.I)  -  L41  •  C(1 .4.  I)  -  L42  •  C(2.4.I) 


C(1  .5.1)  -  C(1 ,5.1)  •  L1I 


-  L43  •  C(3.4.I))  .  L4I 
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6 
1000 


C 

c 
c 


C(2.5.I) 
C(3.5.I) 


(C(2.5.I 
(C(3.5.I 


C(4.5.I)  -  (C(4.5.I)  ■ 

C(3.1 ,1)  - 
C(2.1 .  I)  - 

C(1  .1  .1)  -  C(1 ,1  ,1)  - 


L21  •  C(1 .5, 

L31  •  C(1 ,5. 

L32  •  C(2.5. 

L41  •  C(1 .5. 


C(3.1.I) 
C(2.1.I) 


))  •  L21 
) 

))  •  L3I 

)  -  L42  •  C(2.5.I) 
-  L43  •  C(3.5. I))  .  L4I 


C(3,2.I) 
C(2.2.I) 

C(1  .2.1) 

C(3.3. I) 
C(2.3.I) 


C(3.2.I)  - 

C(2.2. I)  - 

C(1 .2,1)  - 

C(3.3, I)  - 

C(2.3.I)  - 


C(1 ,3.1)  -  C(1 .3,1)  - 


C(3.4,I) 
C(2.4,I) 


C(3.4.I)  - 
C(2.4,  I)  - 


C(1.4.I)  -  C(1 ,4.1)  - 


C(3.5.I) 
C(2,5.I) 

C(1,5.I) 

3  CONTINUE 


C(3.5,I)  - 
C(2,5.I)  - 

C(1 .5.1)  - 


U34 

»   C( 

U24 

»  C( 

U23 

•  c( 

U14 

»  C( 

U13 

►  c( 

U34 

»  C( 

U24 

•  c( 

U23 

»  C( 

U14  ■ 

»  C( 

U13  - 

»  c( 

U34  ■ 

»  C( 

U24  . 

*   C( 

U23  ■ 

*   C( 

U14  « 

>  C( 

U13  < 

>  C( 

U34  . 

>  C( 

U24  « 

'  c( 

U23  < 

.  C( 

U14  < 

-  C( 

U13  < 

'  c( 

U34  « 

>  C( 

U24  « 

1  c( 

U23  < 

'  C( 

U14  < 

'  C( 

U13  < 

'  c( 

3.5.1 


-  U12  •  C(2.1,I) 


-  U12  •  C(2,2,I) 


-  U12  •  C(2.3.I) 


-  U12  •  C(2.4,I) 


-  U12  •  C(2.5.I) 


DO  6 
DO  9 
GG(I1 
DO  6 
DO  6 
HH(I1 
CONTINUE 
CONTINUE 


11 
I 

I 
12 
I 

12 


-1.4 
=  2  ,  IMAX  -  1 
K)     =  C( 1 1 .1.1) 
»  1  .  4 
-  2  .  IMAX  -  1 
I.K)  -  C(I1 ,12+1.1) 


BACKWARD  SUBSTITUTION 


DO  7 
DO  7 
DO  7 
GG(I1 


K 
11 
I 
.1 


«) 


KMAX  -  3, 

1  .  4 

2  .  IMAX  - 
GG(H.I.K) 


1 
2 
3 
7  CONTINUE 

RETURN 

END 


-  1 


f*(I1 
*(I1 
*(I1 
*(I1 


•  K) 
K) 
K) 


I.K) 


GG(1 .I.K+1) 
GG(2.I,K+1) 

GG(3.I  .K-j-1) 
GG(4, I.K+1) 


C 

c- 
c 


C 

c- 
c 


SUBROUTINE  METRIC 

COMMON/F I X/OMEGA 

COMMON/DGRID/DT, IMAX , KMAX. ITEL.  ITEU 

COMMON/GRID1/X(161 , 41 ) , Z( 161 , 41 ) 

COMMON/GR I D/YACOB (161.41) 

COMMON/MTRIX/XIX(161 .41) ,XIZ(161 . 41 ) , ZETAX( 1 61 .41). 
+ZETAZ( 161 .41 ). 
1XIT(161 ,41).ZETAT(161 ,41) 


••   SUBROUTINE  METRIC  COMPUTES  THE  METRICS 
THE  UNSTEADY  COEFFICIENTS  ETAT .  ETC. 


IN  BOTH  DIRECTIONS  AND 
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do  1000  k  -  i 
do  1000  i  -  1  , 

XTAU  -  OMEGA  • 
YTAU  =  OMEGA  • 
PRESENT  SET  UP 


KMAX 

I  MAX 

Z(I.K) 

(-X(I.K)) 
IS  FOR  FLOW 


PAST  AN   AIRFOIL. 


!  !  !  !  ICENTRAL  DIFFERENCES  AT  INTERIOR  POINTS,  TWO-POINT  ONE-i.:ED 
! ! ! ! IDIFFERENCES  DOWNSTREAM.  THREE-POINT  AT  OTHER  OUTER  BOUNDARIES 

IF(I  EQ. 1 .OR. I .EQ. IMAX)  GO  TO  10 

XXI  *  .5  •  (X(I+1 ,K)-X(I-1 ,K)) 

ZXI  -=  .5  •  (Z(I+1 .K)-Z(I-1 .«)) 

GO  TO  15 
10  IF(I  EQ. IMAX)  GO  TO  16 


(X(2.K)  -  X(1 
(Z(2,K)  -  Z(1 


K)) 

K)) 


16 


15 


C! 


17 


18 

20 

i  i  i 


XXI  =  1  .0 

zxi  =  1 .0 

GO  TO  15 
XXI  =  1 .0 
ZXI  -  1 .0 
CONTINUE 
IF(K.EQ.1. 
XZET  =  .5 
ZZET  -  .5 
GO  TO  20 

IF(K.EO.KMAX)  GO  TO  1 8 
XZET  ■=  2.  •  X(I.2)-1  .5 
ZZET  =  2.  •  Z(I,2)  -  1.5  • 
GO  TO  20 

XZET  -  1.5  •  X(I ,KMAX)-2.. 
ZZET  -  1.5  •  Z(I ,KMAX)-2.« 
CONTINUE 

! COMPUTE  JACOBIAN 
YACOBI  =  XXI  •  ZZET  -  XZET  • 
YACOB(I .K)  =  1 .  /  YACOBI 
XIX(I.K)  =  ZZET  .  YACOB(I.K) 
XIZ(I.K)  =  -XZET  •  YACOB(I.K) 
XTAU  «  OMEGA  •  Z(I ,K) 
YTAU  »  -  OMEGA  •  X(I ,K) 


*  (X(IMAX.K)  -  X(IMAX-1 

•  (Z(IMAX.K)  -  Z(IMAX-1 

OR.K.EQ.KMAX)  GO  TO  17 
•(X(I .K+1 )-X(I ,K-1 )) 
•(Z(I ,K+1)-Z(I.K-1)) 


X(I.1)  - 
•  Z(I.1) 


K)) 


-  .5 


X(1.3) 

•  Z(1.3) 


X(I ,KMAX-1)+. 
Z(I ,KMAX-1)+. 


ZXI 


5.X(I 
5«Z(I 


KMAX-2) 
KMAX-2) 


XIT(I .K)  - 
ZETAX(I.K) 
ZETAZ(I ,K) 
ZETAT(I ,K) 
1000  CONTINUE 
RETURN 
END 


XIX(I.K)  •  XTAU  -  XIZ(I.K)  •  YTAU 

-ZXI  »  YACOB(I .K) 

XXI  •  YACOB(I ,K) 

-  ZETAX(I.K)  •  XTAU  -  ZETAZ(I.K)  •  YTAU 


C 

C< 

C 


c 
c 
c 
c 


c 

c 


SUBROUTINE  DISS IP 

COMMON/FLOW/Q1(161 , 41 ) ,02( 1 61 . 41 ) . Q3(  1  61  . 41  ) 

COMMON/PERTR/DQ1(161 . 41 ) ,DQ2( 1 61 , 41 ) , D03( 1 61 

COMMON/DGR I D/DT . IMAX . KMAX , I  TEL . I TEU 

COMMON/GR I D/ Y ACOB (161.41) 

COMMON/DAMP/WW . WW I 

DIMENSION  P(161).EPS(161),DIS1(161 ,4).DIS2(161 


04(161 .41) 
41),D04(161 .41) 


.♦) 


THIS  SUBROUTINE  ADDS  THE  FOURTH  ORDER  DISSIPATION  TERMS  TO  THE 
RIGHT  HAND  SIDE 


IM1 
KM1 
IM2 
KM2 


IMAX 
KMAX 
IMAX 
KMAX 


DO  10  K> 
COMPUTE 
DO    1     I    ■■ 


=2  .  KM1 
SWTICHING 
«  1  .  IMAX 


FUNCTION  BASED  ON  SECOND  DERIVATIVE  OF  PRESSURE 
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P(I) 

DO  2 
IP2  - 
IF(I. 


-  .4  • 

I  -1  . 

1+2 

EQ. IM1) 


(04(1 
I  MAX 


,K)-(Q2(I.K)..2+03(I.K)..2)/(2..Q1(I.K))) 


IP2  -  IMAX 


IM2  -  I  -  2 
IF(I.E0.2)  IM2 
IP1  -  I  +  1 


1 


IF(1 

IM  - 
IF(I 
IF(  I 
IF(I 


IP1 


IMAX 


EQ. IMAX) 

I  -  1 

EQ.1)  IM  -  1 

EQ.1)  IM2  -  1 

EQ. IMAX)  IP2  -  IMAX 
EPS(I)  -  ABS(P(IP1)-2..P(I)+P(IM))/ABS(P(IP1)+2.«P(I)+P(IM)) 
VOL  -  2.   /(YACOB(I . K)+YACOB( IP1 .K)) 
VOL  =  1 . 

DIS1 (1.1)  -   (01(IP1 .K)-Q1(I.K)).V0L 
DIS1 (1.2)  -   (Q2(IP1 .K)-Q2(I.K)).V0L 
DIS1 (1.3)  -   (Q3(IP1 .K)-Q3(I.K)).V0L 
HP1  =  Q4(IP1 ,K)+P(IP1) 
HP  -  04(1 ,K)+P(I) 
HM1  =  Q4(IM,K)  +  P(IM) 
HP2  ■=  Q4(IP2.K)  +  P(IP2) 
HPM  =  Q4(IM.K)+P(IM) 


DIS1 ( I ,4) 
DIS2(I .1) 
DIS2(I.2) 

DIS2(I.3) 
DIS2(I .4) 
2  CONTINUE 
DO  15  I  - 

D2P  = 
C22  = 
C11  = 
C22  «= 
C1 1  - 
C! ! ! ! (SWITCH 


(HP1-HP)»V0L 
(01(IP2.K)-3. 
(Q2(IP2.K)-3. 
(03(IP2.K)-3. 


•(Q1(IP1 
•(Q2(IP1 
•(Q3(IP1 


(HP2-3.»(HP1-HP)-HPM)»VOL 


K)-Q1  ( I .K))-Q1 ( IM.K) ).VOL 
K)-Q2(I .K))-Q2(IM.K)).V0L 
K)-Q3(I ,K))-Q3( IM.K)). VOL 


1 


IM1 

AMAX1(EPS(I).EPS(I+1)) 
66.  •  D2P 

AMAX 1(0.0. (1 .-C22)) 
C22  •  WW/YAC08(I  .*) 
C1 1  •  WW/YACOB(I ,K) 
ON  SECOND-ORDER  AND 


DT 
DT 
SWITCH 


C!  !  !  ! 


15 


! IN  VICINITY 
DISI(I.I)  = 
DIS1 (I .2)  - 

DISKI.3)  = 
DIS1(I.4)  -= 
CONTINUE 
DO  16  I  =  2 
DQ1(I.K) 
DQ2(I.K) 
DQ3(I.K) 
D04(I.K) 
CONTINUE 
CONTINUE 
K  DIRECTION 
C! ! ! ! ! FOURTH-ORDER 
DO  30  I  »  2 
WT=  0.5  •  DT 


OF  SHOCKS 


C11  • 

C11  • 

C11  • 

C11  • 

.  IM1 
0IS1 (  I 
DIS1 (  I 
DIS1 ( I 


DIS2(I.1) 
DIS2( I .2) 
DIS2(I.3) 
DIS2(I.4) 


C22 
C22 
C22 
C22 


OFF  FOURTH-ORDER  DISSIPATION 

•  DIS1 (1.1) 

•  DIS1(I.2) 

•  D1S1 ( I .3) 

•  DIS1 (1.4) 


1) 
2) 
3) 


16 
10 


DIS1 ( I .4)  - 


DIS1 ( 1—1 .1) 
DIS1 ( 1—1 .2) 
DIS1 ( 1—1 
DIS1 ( 1-1 


11 


DISSIPATION  ONLY 

IM1 
.  WW  /  YAC0B(I.2) 


W3  =  0.5  •  DT 

DQ1(I .2)  -WT. 
1+DQ1(I.2) 

DQ2(I .2)  -WT. 
1+DQ2(I.2) 

003(1.2)  -WT. 
1+DQ3(I .2) 

DQ4(I,2)  -WT. 
1+D04(I .2) 

WT=  W3 

DQ1 ( I .KM1 )  -WT« 
1+DQ1 (I ,KM1) 

DQ2( I ,KM1)  =WT. 
1+DQ2(I .KM1) 


WW  / 


(01(1.1) 
(02(1.1) 
(03(1.1) 
(04(1.0 


YACOB(I ,KM1) 
2. 


01(1.2)  +  01(1.3)) 
02(1.2)  +  02(1.3)) 
03(1.2)  +  Q3(1.3)) 
-  2.  •  04(1.2)  +  04(1.3)) 


-  2 


-  2 


(Q1(I.KM2) 
(02(1 ,KM2) 


-  2. 


-  2. 


01(1 .KM1)  +01(1 .KMAX)) 
Q2(I.KM1)  +  Q2(I .KMAX)) 
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D03(I.KM1)  -ITT*  (03(I.KM2)  -  2.  • 
1+003(1 ,KM1 ) 

D04(I.KM1)  -WT»  (04(1. KM2)  -  2.  • 
1+0Q4(I.KM1) 

DO  35  K  -  3  .  KM2 
WT-  -  DT  •  WW  /  YACOB(I.K) 
D01(I.K)  -WT.  (01(1 .K+2)  ■ 
1  I.K)  -  4.  •  01 ( I .K-1)  + 
DQ2(I.K)  -WT. (02(1. K+2) 

•  02(1 .K-1)  + 
WT.(Q3(I ,K+2) 

•  03(1 ,K-1)  4 
WT.(04(I .K+2) 
»  04(1 .K-1)  + 


Q3(I.KM1) 
04(1. KM1) 


+  03(1 .KMAX)) 
+  04(1 .KMAX)) 


1  I .K)  -  4 

DQ3(I ,K) 

1 I.K)  -  4 

DQ4(I ,K) 

1  I  .K)  -  4 

35  CONTINUE 

30  CONTINUE 


4.  •  01 (I .K+1 )  +  6.  •  01 ( 
01(1 .K-2))+DQ1(I.K) 

-  4.  •  02(1 .K+1)  +  6.  •  02 ( 
02(1 .K-2))+DQ2(I.K) 

-  4.  .  03(1 .K+1)  +  6.  •  Q3( 
Q3(I.K-2))+DQ3(I ,K) 

-  4.  .  04( I .K+1 )  +  6.  •  Q4( 
04(1 ,K-2))+DQ4(I,K) 


RETURN 
END 


.41). 


41) 


SUBROUTINE  WALLBC 

COMMON/SURF/PSUR (161) 

C0MM0N/GRID1/X(161  ,41).Z(161.41) 

COMMON/PAR/GAMrfA  .  REYREF  .  ALFA  .  ALFA1  .  REDFRE .  AMI  NF .  ALFAI 

COMMON/DGR I D/DT . IMAX.KMAX. ITEL. ITEU 

C0MM0N/GRID/YAC0B(161 ,41) 

COMMON /DAMP/WW , WW I 

C0MM0N/MTRIX/XIX(161  . 41  ) . XIZ( 161 . 41 ) .ZETAX(161 
+ZETAZ(161 .41). 
1XIT(161 .41 ).ZETAT(161 .41) 

C0MMON/FLOW/Q1(161  ,41 ) ,02(161 ,41 ) ,Q3( 161 ,41 ) .04(161 

DIMENSION  C1(4) 

DIMENSION  A(2.2).RHS(2) 
C! ! ! ! IAPPLY  BOUNDARY  CONDITIONS  ON  THE  CUT  AND  THE  AIRFOIL  SURFACE 

DO  9  I-ITEU.IIIAX 

11  =  IMAX  +1-1 

01(1.1)  -  .5  •  (01(1.2)401(11.2)) 

02(1.1)  -.5*(Q2(I.2)+Q2(I1.2)) 

03(1,1)  -  .5  •  (03(1.2)  +  03(11.2)) 

04(1.1)  =  .5  •  (04(1,2)404(11.2)) 

01(11 , 1)=Q1(I.1) 

02(11 . 1 )=Q2( 1.1) 

03(11 .1)=Q3(I.1) 

04(11 ,1)=Q4(I.1) 
9  CONTINUE 

DO  1  1=   ITEL  .  ITEU 

K  =  3 

C1(1) 

C1(2) 

C1(3) 

UC0N3 
1/QUI.K) 

K  -  2 

C1(1) 

C1(2) 

C1(3) 

UC0N2 
1/01(1. K) 

RHS(1)  -  2.  •  UC0N2  -  UC0N3  -  XIT(I,1) 
C      FOR  VISCOUS  FLOWS  SET  UCON  TO  ZERO  ALSO 

IF(REYREF.GT.e.)  RHS(1)  -  -  XIT(I,1) 

A(1  ,1)  -  XIX(I.I) 

A(1.2)  =  XIZ(I.I) 

A(2.1)  =  ZETAX(IJ) 

A(2.2)  =  ZETAZ(I.I) 


XIT(I.K) 
XIX(I.K) 
XIZ(I.K) 
(Q2(I.K).C1(2)+03(1.K).C1(3)) 


XIT(I.K) 
XIX(I.K) 
XIZ(I.K) 
(02(I.K).C1(2)+03(I.K).C1(3)) 
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c 

c« 
c 


RHS(2)  -  -  ZETAT(I .1) 
TEMPI  -  A(1.1) 
TEMP2  -  A(1.2) 
TEMP3  -  A(2.1) 
TEMP*  -  A(2.2) 

DEN  -  1.  /(TEMPI  •  TEMP4  -  TEMP2  •  TEMP3) 
A(1  .1)  -  A(2.2)  •  DEN 
A(1  .2)  -  -  TEMP2  •  DEN 
A(2, 1)  -  -  TEMP3  •  DEN 
A(2.2)  -  TEMPI  •  DEN 
01(1  .1)  -  2.  •  01(1 .2)  -  01(1,3) 
02(1.1)  ■=  Q1(I.1)*(A(1,1)*RHS(1)+A(1.2)*RHS(2)) 
03(1.1)  «  01(1.1)  •(A(2.1)«RHS(1)+A(2.2).RHS(2)) 
1  CONTINUE 

DO  19  I=ITEL  .ITEU 
U2=02(I.2)/01(I.2) 
W2=O3(I,2)/01(I .2) 

P2=(GAMMA-1 . )«(Q4(I ,2)-0.5»Q1(I . 2)* (U2»U2+W2«W2) ) 
U3=02(I.3)/01(I,3) 
W3=Q3(I.3)/Q1(I,3) 

P3=(GAMHA-1 . )»(Q4(I ,3)-0.5*Q1(l .3)»(U3»U3+W3«W3)) 
P1=(4. .P2-P3)/3. 

PSUR(I)«(GAMMA»P1-1 . )/(.7»AMINF»»2) 
U1=02(1.1)/01(I.1) 
W1=03(I ,1)/01(I .1) 
10  04(1 . 1)«P1/(GAMMA-1 . )+«.5»Q1(I . 1 )• (U1 »U1+W1 .W1 ) 
RETURN 
END 


SUBROUTINE  STRESS( ITN.DALFA) 

C0MM0N/FL0W/Q1(161 . 41 ) .Q2( 1 61 , 41 ) ,Q3( ' 61 .41). 04(1 61 ,41) 

COMMON/DGRID/DT. IMAX.KMAX. I  TEL. ITEU 

C0MM0N/GRID1/X(161 ,41).Z(161.41) 

COMMON/PAR/GAMMA. REYREF. ALFA, ALFA1 .REDFRE . AMINF . ALFAI 

C0MM0N/PERTR/DQ1 (161 ,41),DQ2(161 .41) .003(161 .41 ) .DQ4( 161 .41 ) 

C0MM0N/MUTUR/CMU(161  .41) 

DIMENSION  AA(161 .41) 
1 .RH1(161 ) ,RH2(161) ,RH3(161 ) ,RH4(161) 

COMMON/LOG IC/RSTRT . P I TCH , PLUNGE 

LOGICAL  RSTRT, PI TCH. PLUNGE 

U(I.J)  =  02(1. J)  /  01(1. J) 

V(I .J)  =  03(1. J)  /  01(1. J) 
C     THIS  SUBROUTINE  ADDS  VISCOUS  TERMS  TO  THE  RIGHT  HAND  SIDE 

GOGM  =  GAMMA  /  (GAWMA  -  1 . ) 

IF(ITN.GT. 10. OR. (RSTRT))  CALL  EDDY(DALFA) 
C      COMPUTE  U  AND  V 

KMAXM1  =  KMAX  -  1 

IMAXM1  =  IMAX  -  1 

PR  =  1  . 

DO  10  K  -  1  ,  KMAX 

DO  10  I  -  1  .  IMAX 

E      -  04(1. K)  /  QI(I.K)  -  0.5  •  (U(I ,K)».2+V(I.K).*2) 
10  AA(I ,K)  »  GOGM  .  E 
C 

C     COMPUTE  TXX.TXY  AND  VISCOUS  DISSIPATION  AT  I  -  1  /  2 
C 

DO  30  K  -  2  .  KMAXM1 

DO  20  I  »  2  .  IMAX 

UXI  «  U(I ,K)  -  U(I-1 ,K) 

VXI  =  V(I .K)  -  V( 1-1 ,K) 

AXI  =  AA(I ,K)  -  AA(I-1 .K) 

UZET=  .25* (U( I .K+1 )-U( I .K-1 )+U( 1-1 .K+1 )-U( 1-1 ,K-1 ) ) 

VZET=  .25»(V(I ,K+1)-V(I,K-1)+V(I-1,K+1 )-V( 1-1 .K-1)) 

AZET=  ,25«(AA(I ,K+1)-AA(I ,K-1)+AA(I-1 ,K+1)-AA(I-1 .K-1)) 

XXI  «  X(I ,K)  -  X( 1—1 ,K) 
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ZXI  «=  Z(I.K)  -  Z(I-1  .K) 

XZET-  .25  •  (X(l ,K+1)-X(I ,K-1)+X(I-1 ,K+1 )-X( 1-1 ,K-1 ) ) 

ZZET-  .25  •  (Z(I.K+1)-Z(I.K-1)+Z(I-1.K+1)-Z(I-1.K-1)) 

YAC  -  XXI  •  ZZET  -  ZXI  •  XZET 

YAC  -  1 .  /  YAC 

XIX  -  ZZET  .  YAC 

ZETAX-  -  ZXI  •  YAC 

XIZ  -  -XZET  »  YAC 

ZETAZ-  XXI  •  YAC 

CNM  =  .5  •  (CMU(I.K)  +  CMU(I-I.K)) 

UX   =  UXI  »  XIX  +  UZET  •  ZETAX 

VX   -  VXI  •  XIX  +  VZET  •  ZETAX 

AX   «=  AX  I  •  XIX  +  AZET  •  ZETAX 

UZ   =  UXI  •  XIZ  +  UZET  •  ZETAZ 

VZ   =  VXI  •  XIZ  +  VZET  «  ZETAZ 

AZ   -  AXI  •  XIZ  +  AZET  •  ZETAZ 

TXX  -  -(-•*.  •  UX  +  2.  •  VZ)  •  CNM  /  3. 

TXY  =  CNM  .  (UZ  +  VX) 

TYY  -  -CNM  /  3.  •  (-4.  •  VZ  +  2 .  •  UX) 

R4   =  ((U(I ,K)+U(I-1 ,K))»TXX+(V(1 ,K)+V(I-1 , K)  )»TXY)«0. 5 
1      +  CNM  /  PR/ (GAMMA  -  1 . )  •  AX 

S4   =  ((U(I ,K)+U(I-1 ,K))*TXY+(V(I ,K)+V(I-1  .  K)  )  .TYY) .0 . 5 
1      +  CNM  /  PR  /  ( GAMMA  -  1  .  )  •  AZ 
C       DEBUG 

C       TURN  OFF  ENRGY  DISSIPATION  AND  DIFFUSION 
R4  =  0. 
S4  =  0. 
RH1(I)  =  0. 

RH2(I)  =  (XIX  •  TXX  +  XIZ  •  TXY)  /  YAC 
RH3(I)  -  (XIX  »  TXY  +  XIZ  •  TYY)  /  YAC 
20  RH4(I)  -  (XIX  •  R4  +  XIZ  •  S4)  /  YAC 
DO  30  I  =  2  .  IMAXM1 

DQ1(I,K)  -  D01(I,K)  +  RH1(I+1)  -  RH1(I) 
DQ2(I.K)  -  DQ2(I.K)  +  RH2(I+1)  -  RH2(I) 
DQ3(I,K)  -  DQ3(I.K)  +  RH3(I+1)  -  RH3(I) 
D04(I,K)  «=  D04(I.K)  +  RH4(I+1)  -  RH4(I) 
30  CONTINUE 
C      IN  THE  Z  DIRECTION 

DO  70  I  =  2  ,  IMAXM1 

DO  60  K  =  2  ,  KMAX 

UXI  =  .25  •  (U(  I-t-1  ,K)-U(  1-1  .K)+U(  1  +  1  ,K-1  )-U(  1-1  .K-1  )) 

VXI  =  .25  •  (V(l  +  1 ,K)-V(I-1 ,K)+V(I  +  1 ,K-1)-V(I-1  .K-1)) 

AXI  =  .25  •  (AA(I  +  1 ,K)-AA(I-1  ,K)+AA(I  +  1 ,K-1)-AA(I-1  .K-1)) 

XXI  »=  .25  •  (X(I+1 ,K)-X(I-1 .K)+X(I+1 ,K-1)-X(I~1.K-1)) 

ZXI  =  .25  •  (Z(I  +  1  ,K)-Z( 1-1 . K)+Z( 1  +  1 .K-1 )-Z( 1-1 .K-1  )) 

UZET  =  U(I  ,K)  -  U(I ,K-1) 

VZET  =  V(I,K)  -  V(I ,K-1) 

AZET  =  AA( I .K)  -  AA(I ,K-1) 

XZET  =  X(I  ,K)  -  X(I  .K-1) 

ZZET  =  Z(I ,K)  -  Z(I ,K-1) 

YAC  =  XXI  •  ZZET  -  ZXI  »  XZET 

YAC  =  1 .  /  YAC 

XIX  =  ZZET  •  YAC 

ZETAX-  -  ZXI  •  YAC 

XIZ  -  -XZET  •  YAC 

ZETAZ=  XXI  •  YAC 

CNM  =  .5  •  (CMU(I.K)  +  CMU(I.K-I)) 

UX   -  UXI  .  XIX  +  UZET  »  ZETAX 

VX   =  VXI  •  XIX  +  VZET  •  ZETAX 

AX   =  AXI  •  XIX  +  AZET  •  ZETAX 

UZ   =  UXI  »  XIZ  +  UZET  •  ZETAZ 

VZ  =  VXI  •  XIZ  +  VZET  •  ZETAZ 

AZ   =  AXI  •  XIZ  +  AZET  •  ZETAZ 

TXX  =  -(-4.  •  UX  +  2.  •  VZ)  •  CNM  /  3. 

TXY  =  CNM  •  (UZ  +  VX) 

TYY  =  -CNM  /  3.  •  (-4.  •  VZ  +  2 .  •  UX) 

R4   =  ((U(I ,K)+U(I ,K-1)).TXX+(V(I ,K)+V(I ,K-1)).TXY).0.5 
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1      +  CNM  /  PR/(GAMMA  -  1  .  )  •  AX 

S4   -  ((U(I.K)+U(I.K-1)).TXY+(V(I.K)+V(I,K-1)).TYY)»0.5 
1      +  CNM  /  PR  /  (GAMMA  -  1 . )  •  AZ 

R4  -  0. 

S4  -  0. 

RH1 (K)  -  0. 

RH2(K)  -  (ZETAX  •  TXX  +  ZETAZ  •  TXY)  /  YAC 

RH3(K)  =  (ZETAX  •  TXY  +  ZETAZ  •  TYY)  /  YAC 
60  RH4(K)  -  (ZETAX  •  R4  +  ZETAZ  •  S4)  /  YAC 

DO  70  K  «  2  .  KMAXM1 

D01(I. K)  =  D01(I.K)  +  RH1(K+1)  -  RH1(K} 

DQ2(I.K)  =  DQ2(I,K)  +  RH2(K+1)  -  RH2(K) 

DQ3(1.K)  ■  DQ3(I,K)  +  RH3(K+1)  -  RH3(K) 

D04(1.K)  -  DQ4(I.K)  +  RH4(K+1)  -  RH4(K) 
70  CONTINUE 

RETURN 
END 


C 

c< 

c 


SUBROUT I NE  LOAD (CPS . CL . CD . CM , ALFAS) 

C0Mw(0N/GRID1/X(161  .41).Y(161.41) 

COMMON/DGR I D/DT . IMAX.KMAX. ITEL. ITEU 

DIMENSION  CPS(161) 
C 

C      THIS  SUBROUTINE  COMPUTES  THE  INVISCID  CONTRIBUTIONS 
C     TO  LOADS  ON  THE  AIRFOIL  SURFACE 
C 

IMAXM1  -  IMAX  -  1 

CL  =  0. 

CD  -  0. 

CM  =  0. 

DO  400  I  -  ITEL  .  ITEU  -  1 

XL  =  .5  •  (X(I.1)+X(I+1 ,1)) 

YL  -  .5  •  (Y(I,1)+Y(I+1  .1)) 

DX  -  X(I  +  1  .1)  -  X(I,1) 

DY  -  Y(I+1 .1)  -  Y(I,1) 

CPA  «  CPS(I+1)  •  .5  +  CPS(I)  •  .5 

DCL  =  CPA  •  (-DX) 

DCD  -  CPA  •  DY 

CL  -  CL  +  DCL 

CD  =  CD  +  DCD 
400  CM  =  CM  +  DCD  •  YL  -  DCL  •  XL 
C 

DCL  -  CL  •  COS(ALFAS)  -  CD  •  SIN(ALFAS) 

CD  =  CL  •  SIN(ALFAS)  +  CD  •  COS(ALFAS) 

CL  =  DCL 

RETURN 

END 


C 

c« 

c 


SUBROUT I NE  WRAP ( 1 1 . J  J . XS I NG . YS I NG . XP , YP , S0 . A0 . Y0 ) 
DIMENSION  S0(161 ,4).YB(41 ,4).A0(161 .4) ,XP( 1 ) . YP( 1 ) 

C 

C      THIS  SUBROUTINE  UNWRAPS  THE  AIRFOIL  AND  STORES  THE  UNWRAPPED 

C      SURFACE  IN  ARRAYS  A0  AND  S0 .  IT  ALSO  DETERMINES  THE  STRETCHING 

C      IN  THE  ETA  DIRECTION. 


C 


IMID  -  (II  +  1)  /  2 
DY   -  .8  /  (JJ  -  2) 
DO  1  J  =  2  .  JJ 
Y  =  FL0AT(J-2)  •  DY 
1  Y0(J.1)  =  1 .25  •  Y  /  (1 .  -  Y  •  Y) 
Y0(1  ,  1)  -  -  Y0(3.1) 
PI  .  4.  •  ATAN  (  1 . ) 
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ANGL  -  PI  +  PI 

U  m   XP(1)  -  XSING 

V  -  YP(1)  -  YSING 
U  -  1. 

V  »  0. 

I IM1  -  II  -  1 

DO  2  I  -  1  .  II 

X11  «=  XP(I)  -  XSING 

Y11  =  YP(I)  -  YSING 

ANGL  «=  ANGL  +  ATAN2(  (U*Y1 1-V.X1  1  )  .  (I>X1 1+V.Y1 1  )  ) 

R   -  SQRT(X11«»2  +  Y11»«2) 

U    =  X1  1 

V  =Y11 

R  =  SORT(R) 

A0(I.1)  =  R  •  COS(.5  •  ANGL) 
2  S0(I.1)  -  R  •  SIN(.5  •  ANGL) 
C! ! ! ! ! IF  OUTPUT  OF  UNWRAPPED  COORDINATES  IS  DESIRED 
C       WRITE  (6.1000) 

C      WRITE  (6.2000)  (I ,A0(I . 1) ,S0(I  .  1) . I  -  1  .  II) 
RETURN 
1000  FORMAT (1X. 'UNWRAPPED  COORDINATES  IN  THE  TRANSFORMED  PLANE') 
2000  FORMAT(I5  ,  2F20.8) 
END 


SUBROUTINE  TABINT(XP, YP. XSING. YSING, N) 
DIMENSION  XP(161) ,YP(161 ) ,S0(161 ) ,A0(161) 
C! ! ! ! ! SMOOTH  THE  AIRFOIL  SURFACE  BY  FINDING  ADDITIONAL  POINTS 
U  =  XP(1)  -  XSING 

V  =  YP(1)  -  YSING 
U  =  1  . 

V  -  0. 

ANGL  -  8.  •  ATAN(1 .) 

DO  1  I  =  1 . N 

X11  =  XP(I)  -  XSING 

Y11  =  YP(I)  -  YSING 

ANGL  =  ANGL  +  ATAN2 ( (U* Y1 1-V»X1 1 ) . (U«X11+V»Y1 1 ) ) 

R  =  SORT(X1 1**2  +  Y11  ••  2) 

U  =  X11 

V  =  Y11 

R  =  SORT(R) 

A0(I)  =  R  •  COS(ANGL  •  .5) 
1  S0(I)  =  R  •  SIN(ANGL  •  .5) 
DX  =(A0(N)-A0(1))/96. 
A0ST  =  A0(1 ) 
DO  3  I  =  1  .97 
XX  =  FLOAT(I-I)  .  DX  +  A0ST 
CALL  TAINT(A0.S0,XX.YY.N.3.NER.MON) 
XP(I)  =  XX  •  XX  -  YY  .  YY  +  XSING 
3   YP(I)  -  2.  •  XX  •  YY  +  YSING 
RETURN 
END 


SUBROUT I NE  TA I  NT ( XTAB . FTAB . X . FX , N , K . NER . MON ) 
DIMENSION  XTAB(1) .FTAB(1) ,T(10),C(10) 

C 

C      NASA  -  AMES  SUBROUTINE  FOR  POLYNOMIAL  INTERPOLATION 

C     OF  A  TABULATED  FUNCTION 


C 


IF(N-K)  1  ,  1 

1  NER  =  2 
RETURN 

2  IF(K-9)  3,3.1 

3  IF(MON)  4,4,5 
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5  IF(M0N-2)  6.7,4 
4  J  -  0 

NM1  ■  N  -  1 

DO  8  I  -  1  ,  NM1 

IF(XTAB(I)  -  XTA8(I+1))  9.11.10 

1 1  NER  -  3 
RETURN 

9  J  -  J-1 

GO  TO  8 
10  J  =  J+1 
8  CONTINUE 

MON  =  1 

IF(J)  12  .  6  ,  6 

12  MON  =  2 

7  DO  13  I  -  1  .  N 

IF(X  -  XTAB(I))  14.14.13 

14  J  =  I 

GO  TO  18 

13  CONTINUE 
GO  TO  15 

6  DO  16  I  -  1  .  N 
IF(X-XTAB(I))  16.17.17 

17  J  -  I 

GO  TO  18 
16  CONTINUE 

15  J  -  N 

18  J  -  J  -  (K+1)  /  2 
IF(J)  19.19.20 

19  J  ■  1 

20  M  -  J  +  K 

IF(M  -  N)  21 .21 .22 

22  J  =  J  -  1 
GO  TO  20 

21  KP1  -  K  +  1 
JSAVE  «=  J 

26  DO  23  L  -  1 .  KP1 
C(L)  =  X  -  XTAB(J) 
T(L)  -  FTAB(J) 

23  J  -  J+1 

DO  24  J  -  1 ,K 
I  «=  J+1 
25  T(I)  =  (C(J).T(I)-C(I).T(J))/(C(J)-C(I)) 
I  =  1  +  1 
IF(I-KP1)  25.25.24 

24  CONTINUE 

FX  -  T(KP1) 
NER  -  1 
RETURN 
END 


SUBROUTINE  SING(N2.N.X.Z.XLE. YLE.TEA.TES .XSING. YSING.CHD) 
C 
C 
C      THIS  SUBROUTINE  COMPUTES  SINGULAR  POINT  LOCATIONS. 


DIMENSION  X(2)  .  Z(2) 

NU  =  N2 

N1  -  N2  +  1 

N3   -  N2  -  1 

X1   =  X(N1) 

Z1   =  Z(N1) 

X2   =  X(N2) 

Z2   -  Z(N2) 

X3   =  X(N3) 

Z3   -  Z(N3) 

D1   -  X2  ••  2  -  X1  •• 
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D2 

-  Z2 

•  •  2 

-  Z1 

•  • 

2 

D3 

-  X2 

-  XI 

D4 

-  Z2 

-  Z1 

D5 

-  X3 

•  •  2 

-  X1 

•  • 

2 

D6 

-  Z3 

•  •  2 

-  Z1 

•  • 

2 

D7 

-  X3 

-  XI 

D8 

-  Z3 

-  Z1 

B 

-=  (D7 

•  ( 

D1  + 

D2) 

-  D3< 

>(D5+D6))/(2.' 

»(D7« 

■D4-D3' 

►  08)) 

IF(ABS(D3).LT. 

ABS(D7)) 

GO  TO  10 

A  « 

(D1  + 

D2  - 

■  2.  < 

«  B 

•  D4) 

/  (2- 

•  D3) 

GO 

TO  20 

10 

A  = 

(D5  + 

D6  - 

■   2.  « 

'  B 

•  D8) 

/  (  2. 

•  07) 

20 

CONTINUE 

R  = 

SORT( 

(X2-/ 

»)••  : 

!  + 

(Z2-B) 

•  •2) 

XLE  -  X(NU) 

YLE  =  Z(NU) 

CHD  -  X(1)  -  X(NU) 

A2  =  (Z(2)-Z(1))  /(X(2)  -  X(1)) 

A1  -  (Z(N)-Z(N-1))/(X(N)-X(N-1)) 

TES  =  .5  •  (A1  +  A2) 

TEA  =  A2  -  A1 

TEA  =  TEA  •  57.29578 

XSING  =  (A+XLE)  /2. 

YSING  =  (B+YLE)  /  2. 

RETURN 

END 


SUBROUTINE  AIRFOL( 1 1 . JJ . IT . IE) 
C0MM0N/GRID1/X(161  .41 )  ,Z(161  .41 ) 
COMMON /YSYM/ ISYM 

DIMENSION  S0(161 ,4).A0(161 .4).Y0(41 . 4) . XP( 1 61 ) . YP( 1 61 ) . 
1E(161).F(161).B0(49) 
C 

DATA  (B0(I).I«1 ,32)/1 . .1 .0414,1 .0836.1 . 1270. 1 . 1715. 1 . 2175. 1 .2651 . 
11 .3145.1 .3659.1 .4199.1 .4755.1 .5349.1 .5973,1 .6636.1 .7342.1 .8099. 
21 .8914.1 .9799.2.0764.2.1829.2.3012.2.4341 ,2.5653.2.7597.2.9646. 
33 . 2 1 06 . 3 . 51 41 . 3 . 90 1 9 . 4 . 42 1 9 . 5 . 1 687 . 6 . 3632 . 8 . 6809/ 
C 

C! ! ! ! iCOMPUTE  THE  COMPUTATIONAL  GRID  POINTS  BASED  ON  INPUT  AIRFOIL  SHAPE 
DOS  I  *  1  .32 
8  A0(I .1)  -  B0(I) 
READ  (5,1) 

1  FORMAT (1X) 

READ  (5.2)  FNU.FNL.ZSYM 

2  FORMAT(3F10.0) 
ISYM  =  0 
IF(ZSYM.NE.0.)  ISYM  -  1 

•       II  =  157 
«       JJ  =   41 


IT  = 

=   31 

IE  i 

=  127 

IIP1 

=  II 

+  1 

IIM1 

=  II 

-  1 

IIJJ 

-  II 

•  JJ 

IIJJ2  =  ] 

•  (  JJ- 

-2) 

ILE 

-  (IT 

+  IE  ) 

/  2 

ISTP 

=  0 

NN 

=  5 

NRF 

-  0 

NOTAPE  -  1 

PI 

=  4. 

•  ATAN(' 

I.) 

NU  = 

FNU 

NL  - 

FNL 

N   = 

NU  + 

NL  -  1 

READ(5.1) 

READ  (5.333)  (XP( I ) . YP( I ) . 1  «  NL  ,  N) 
333  FORMAT(2F10.0) 

9994  FORMAT(F20.8) 
L  «  N  +  1 

IF(ZSYM  .GT.  0. )  GO  TO  9995 

L  ■  NL  +  1 

READ(5.1) 

READ  (5.333)  (XP(L-I) .YP(L-I) , I«1 ,NL) 

GO  TO  9996 

9995  K1  -  L 

DO  16  1  -  NL  ,  N 
K  =  K1  -  I 
XP(K)  =  XP(  I) 
YP(K)  -  -  YP( I ) 
16  CONTINUE 

9996  SCALE  =  1.  /(XP( 1 )-XP(NL) ) 
XX     ■=  XP(NL) 

YY     =  YP(NL) 
DO  9997  I  =  1  ,  N 
XP( I)  =  XP( I)  •  SCALE 

9997  YP( I)  =  YP(I)  •  SCALE 

CALL  SING(NU,N,XP,YP,XLE.ZLE.TEA.TES,XSING.YSING.CHD) 
CALL  TABINT(XP.YP.XSING,YSING,N) 
NBODY  =  IE  +  1  -  IT 
DO  6791  1=1  ,  NBODY 
L  -  I  -  1 
E(IT+L)  =  XP(I) 
6791  F(IT+L)  =  YP(I) 
IEP1  =  IE  +  1 
SLOPT  m   TES  •  .75 
DO  438  I  =  IEP1  .  II 
11  -  I  +1  -  IE 
E(I)  =  A0(I1  ,1) 
DXI  *  1 .  /  48. 

D    =  4.  /  3.  •  (E(I)  -  .25) 
F(I)  =  F(IE)  +  SLOPT  .  ALOG(D)  /  D 
L  =  IIP1  -  I 
E(L)  =  E(I) 

438  F(L)  =  F(IT)  +  SLOPT  .  ALOG(D)/D 
C      WRITE  (6.439) 

439  FORMAT (2X.3H   I . 19X . 1 HX . 19X . 1HY) 

C      WRITE  (6.37)  ( I . E( I ) . F(  I )  .  1  =  1  ,  II) 

CALL  WRAP  01 . JJ.XSING.YSING.E.F.S0.A0.Y0) 
C! ! ! ! !MAP  GRID  BACK  TO  PHYSICAL  PLANE  AND  SHIFT  TO  QUARTER  CHORD 
DO  10  J  =  2  .  J  J 
DO  1 0  I  =  1  ,  II 

X(I.J-1)  =  A0(I,1)«.2  -  (S0(I .1 )+Y0(J.1))..2 
1  -  0.25 
10  Z(I.J-I)  =  2.  •  A0(I.1)  .  (S0(I  ,1)+Y0(J,1)) 
JJ  =  JJ  -  1 
RETURN 
37  FOPMAT(I5.2F20.8) 
END 


C 


SUBROUTINE  CLUSTR(DS) 

COMMON/GRID1/X(161 ,41).Z(161,41) 

COMMON/DGRID/DT. IMAX.KVUX. I  TEL. ITEU 

DIMENSION  S(41) ,XP(41).YP(41) ,R(41) 
C 

C      THIS  SUBROUTINE  CLUSTERS  A  GIVEN  X.Z  GRID  SUCH  THAT  THE  FIRST  POINT  IS  AT 
C      THE  USER-SPECIFIED  DISTANCE  DNMIN 
C! ! ! ! ICOMPUTE  THE  OLD  STRETCHING 

DO  100  I  =  1  .  IMAX 

SO)  =  0. 

XP(1)  =  X(I.1) 
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YP(1)  -  Z(I.1) 
DO  10  K  -  2  ,  KMAX 
XP(K)  -  X(I ,K) 
YP(K)  -  Z(I ,K) 
10  S(K)  -  SQRT((XP(K)-XP(K-1)).*2+(YP(K)-YP(K-1)).»2) 
1+S(K-1) 
SUMDX  -  S(KMAX) 

CALL  STRTCH ( SUMDX. DS.F1 . KMAX . FACTOR) 
:       WRITE  (6,200)  I .FACTOR 
R(1)  -  0. 
DR  =  DS 

DO  20  K  «  2  ,  KMAX 
R(K)  -  R(K-1 )  +  DR 
DR  =  DR  •  FACTOR 
20  CONTINUE 

RLAST  -  1 .  /  R(KMAX) 
DO  30  K  -  2  .  KMAX 
R1  =  R(«)  •  RLAST  •  S(KMAX) 
!  !  !  !  IREDISTRIBUTE  THE  CONSTANT-ETA  LINES 

CALL  TAINT(S,XP.R1 ,XP1 .KMAX . 3 ,NER .MON) 
X(I  .K)  -  XP1 

CALL  TAINT (S. YP. R1 ,YP1 .KMAX . 3 . NER.MON) 
Z(I  ,K)  -  YP1 
30  CONTINUE 
100  CONTINUE 

WRITE  (6.115) 
DO  110  I  =  1  .  IMAX 
DX  =  X(I .2)  -  X(I 
DY  =  Z(1.2)  -  Z(I 
DN  =  SORT(DX»DX+DY»DY) 
WRITE(6.120)  I  .  DX  ,  DY  .  DN 
110  CONTINUE 

RETURN 
115  FORMAT ( 5X. 6HNORMA L. 1X.8HD I  STANCE. 3H  AT . 4H  THE.5H  WALL./ 

1 . 5H     I . 8X . 2HDX . 8X . 2HDY . 8X . 2HDN . //) 
120  FORMAT(I5.3F10.5) 
200  FORMAT(I5.F10.5) 
END 


!l 


c 

c« 

c 

c 
c 
c 
c 
c 


SUBROUTINE  STRTCH ( SUMDX ,DX1 ,F1 ,N1  ,R) 

THIS  SUBROUTINE  DETERMINES  A  GEOMETRIC 

PROGRESSION  OF  GRID  SPACING  BETWEEN  1  AND  N1  SUCH  THAT 

SUM8DX)  EOUALS  SUMDX.  THE  RATIO  BETWEEN  SUCCESSIVE 

SPACINGS  IS  R. 

N  «  N1  -  1 

R  =  1  .5 

E1 

E2 

DO 

F= 

FP 


=  1 . E-04 
=  1 . E-04 
10  L  =  1 .  50 
(R-1 )  •  SUMDX  - 
=  SUMDX  -  DX1  • 
RITER  =  R  -  F/  FP 
C      IF(1 .E-02.LT. RITER 
C      IF(1  .  .LT. RITER. AND 
IF(ABS(R-RITER).LT 
R  -  RITER 
10  CONTINUE 
R  =  1 .0001 

DX1  =  DZTOT/FLOAT(N1-1) 
RETURN 
1  R=  RITER 
RETURN 
END 


DX1«(R«»N-1 ) 
FLOAT(N)  .  R 


(N-1) 


AND.RITER.LT. 1 .)  RITER  - 
RITER. LT. 100.)  RITER=.01 
R»E1)  GO  TO  1 
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c 


SUBROUTINE  EDDY(DALFA) 

COMMON/F L0W/Q1 (161 .41) .02(161 .41 ).Q3( 161 ,41) ,Q4(161 .41) 

COMMON  /MUTUR/CMU  (161,41) 

COMMON/SK I NCF/CF( 1 6 1 ) 

C0MM0N/DGR1D/DT. IMAX.KMAX. I  TEL. ITEU 

COMMON/PAR/GAMMA. REYREF. ALFA. ALFA 1 . REDFRE . AMI NF . ALFA I 


C0MM0N/GRID1/X(161  ,41),Z(161.41) 

(41) 


DIMENSION  TIN(41 ) ,T0UT(41 ) .Y(- 
C 
C       INITIALIZE  VISCOSITY  EVERYWHERE 

FACT1  =  DT  •  AMINF  /  REYREF 

CMUMAX  =  100.  *  FACT1  /  DT 

DO  1  K  «=  1  .  KMAX 

DO  1  1  =  1  .  I  MAX 
1  CMU( I ,K)  =  FACT1 
C      THIS  SUBROUTINE  COMPUTES  THE  EDDY  VISCOSITY  BASED  ON  THE 
C      BALDW1N-L0MAX  TWO  LAYER  MODEL 
C 

DO  100  I  =  2  .  IMAX  -  1 

UDIF  «  0. 

FMAX  =  0. 1E-06 

YMAX  -  . 1 E-06 

FYMAX  -  0. 

Y(1)  -  0. 

UWALL  =  0. 

IF(I  GT. ITEU.OR. I .LE. ITEL)UWALL  =  SQRT(Q2( I . 1 )»»2+03( I , 1 )..2)/ 
101(1 .1) 
C      COMPUTE  TAU  AT  THE  WALL 

UET  =  1  .(02(I .2)/Q1(I.2)  -  Q2( I . 1 )/Q1 ( I . 1 ) ) 

VET  -  1 ..(Q3(1.2)/Q1(I.2)  -  Q3(I.1)/01(I .1)) 

XXI  =  X(I  +  1 .1  )  -  X(I-1 , 1 ) 

ZXI  =  Z(I+1 .1)  -  Z(I-1 .1) 

XET  =  4.  •  X(1.2)  -  3.  •  X(I.1)  -  X(1.3) 

ZET  =  4.  •  Z( I .2)  -  3.  •  Z(I.1)  -  Z(1.3) 

XXI  =  .5  •  XXI 

ZXI  =  .5  •  ZXI 

XET  =  .5  •  XET 

ZET  =  .5  •  ZET 

YAC  =  1 .  /  (XXI  •  ZET  -  ZXI  •  XET) 

OMEGA  =   (UET  •  XXI  -  VET  .  ZXI  )  •  YAC 

TWALL  =  AMINF  »  OMEGA  /  REYREF 

CF(I)  =  2.  •  TWALL  /  (AMINF. »2) 

FACT  =  S0RT(Q1(I,1)  •  ABS(TWALL)).REYREF/(26.«AMINF) 

DO  10  K  =  2  ,  KMAX-1 

UXI  -  (02(1+1 ,K)/Q1(I+1 ,K)  -  02(1-1 .K)/Q1(I-1 ,K)) 

VXI  =  (Q3(I+1 ,K)/Q1(I+1 ,K)-Q3(I-1 ,K)/01 ( 1-1 .K)) 

UET  =  (02(1 .K+1)/Q1(I.K+1)-02(I.K-1)/Q1(I,K-1)) 

VET  =  (03(1 .K+1)/01(I ,K+1)-Q3(I ,K-1 )/Q1 ( I .K-1 ) ) 

XXI  -  X(I+1 .K)  -  X(I-1 .K) 

ZXI  =  Z(I+1 .K)  -  Z( 1—1 .K) 

XET  =  X(I ,K+1 )  -  X(I ,K-1) 

ZET  =  Z(I.K+1)  -  2(1 .K-1 ) 

YAC  =  1 .  /  (XXI  .  ZET  -  ZXI  •  XET) 

OMEGA  =  ABS(UET»XXI+VET.ZXI-UXI.XET-VXI.ZET)  •  YAC 

UDIF  =  S0RT(Q2(I  ,K). «2+Q3(I  .K)..2)/01(I ,K)  -  UWALL 

IF(ABS(UDIF) .GT.UDIFMAX)  UDIFMAX  =  ABS(UDIF) 

Y(K)  =  SORT((X(I  ,K)-X(I.K-1))..2+(Z(I.K)-Z(I ,K-1))..2)+Y(K-1) 

F  =  Y(K)  •  OMEGA 

IF((Y(K)»FACT) .GT.20.)  GO  TO  31 

I F( I  GT. ITEL.AND. I . LT . ITEU)  F  =  F  .  (1.  -  EXP(-Y(K).FACT) ) 
31  CONTINUE 
C 

C     MODIFIED  TURBULENCE  MODEL  APPLY  FOR  SPECIFIED  RANGE  OF  ANGLES  *HERE 
C      FY  IS  USED  T0  FIND  THE  SECOND  PEAK  VALUE  OF  F  FUNCTION 
C 
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10 


20 


ALFAI .AND. DALFA.GE.0. )  THEN 
Y(K) 

FYMAX)  THEN 
FY 

Y(K) 


IF(ALFA.LT 
FY  -  F  • 
IF(FY.GT 
FYMAX 
FMAX 
YMAX 
END  IF 
END  IF 

I F( ALFA. GE. ALFAI .OR . DALFA. LT . 0 . )  THEN 
IF(F.GT.FMAX)  THEN 
FMAX  ~  F 
YMAX  -  Y(K) 
END  IF 
END  IF 

FCT  -  Y(K)  •  FACT 
IF(FCT.GT.20.  )  FCT  ■=  20. 
FCT  -  ABS(FCT) 
EL  =  .4  •  Y(K)  •  (1.  -  EXP(-FCT)) 


01(1. K)  •  EL 
ABS(TIN(K)) 


TIN(K)  « 

TIN(K)  - 

CONTINUE 

KSWTCH  »  0 

FWAKE  =  YMAX  •  FMAX 

F1  =  0.25  •  YMAX  •  UDIF 

IF(F1 . LT. FWAKE)  FWAKE  « 

DO  20  K  -  2  .  KMAX  -  1 

FKLEB  =  0. 

IF(ABS(Y(K)/YMAX).LT. 1 


EL  •  OMEGA 


•  •2 
F1 


/  FMAX 


FKLEB 
END  IF 
TOUT(K)  - 
TOUT(K)  = 
I F( KSWTCH 
IF(TIN(K) 
CONTINUE 


1 .  /  (1.  +  5.5 


E+04)  THEN 

•  (0.3  •  Y(K)/YMAX)  ••  6) 


.0168  •  1.6  •  Q1(I.K)  •  FWAKE  •  FKLEB 
ABS(TOUT(K)) 
NE.0)  GO  TO  20 
GT.TOUT(K))  KSWTCH  -  K  -  1 


C! ! ! ! ! TOTAL 
DO  30 
IF(K, 


VISCOSITY  IS  SUM 

K  -  2  .  KMAX  -  1 

LE. KSWTCH)  THEN 


OF  LAMINAR  AND  INNER/OUTER  LAYER  AS  APPROPRIATE 


30 
100 


CMU(I .K) 

ELSE 

CMU(I.K) 

END  IF 

CONTINUE 

CONTINUE 

RETURN 

END 


DT 


DT 


(AMINF/REYREF  +  ABS(TIN(K))) 
(AMINF  /  REYREF  +  ABS(TOUT(K) ) ) 


C 
C 
C 

C 


SUBROUTINE  RES I (OMEGA) 

C0MM0N/PERTR/DQ1(161 .41) ,DQ2(161  ,41) .003(161 .41) .004(161 ,41) 

C0MM0N/GRID1/X(161 ,41) ,Z(161 ,41) 

COMMON/DGR I D/DT . I MAX , KMAX . I  TEL , I TEU 

C0MM0N/FL0W/Q1(161 , 41 ) , 02 ( 1 61 , 41 ) ,Q3( 1 61 . 41 ) , 04 ( 1 61 . 41 ) 

COMMON/PAR/GAMMA , REYREF . ALFA . ALFA 1 , REDFRE , AMI NF . ALFAI 

DIMENSION  RHS(161 ,4) 

XTAU(I.K)  «  OMEGA  •  Z(I.K) 

YTAU(I.K)  -  -  OMEGA  •  X(I.K) 

THIS  SUBROUTINE  COMPUTES  THE  RESIDUAL  ON  THE  RIGHT  HAND 

SIDE  ARISING  FROM  THE  EULER-  PART  OF  THE  GOVERNING  EQUATIONS 

FLUX  TERMS  WITHIN  THE  XI-  DERIVATIVE 

DO  100  K  =  2  .  KMAX  -  1 

DO  10  I  =  1  ,  IMAX 

UCON  =  (Q2(I.K)/01(I,K)  )  •  (Z( I ,K+1 )-Z( I ,K-1 ) ) 
1  -  (03(1 ,K)/Q1(I.K))  •  (X(I.K+1)-X(I ,K-1)) 

UCON  =  0.25  •  DT  »  UCON 
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XIT  -  -  XTAU(I.K)  .(Z(I.K+1)-Z(I.K-1)) 
1  +  YTAU(I.K)  •  (X(I.K+1)  -  X( I ,K-1 ) ) 

XIT  -  XIT  .  DT  .  0.25 

UCON  -  UCON  +  XIT 

RHS(I . 1)  -  UCON  •  Q1 (I ,K) 

R  -  1 .  /  01 (I ,K) 

P  =  (GAMMA-1.)  •  (04(1. K)  -  .5  •  R  • (Q2( I ,K) ..2+ 
1  03(1. K).. 2)) 

RHS(I.2)  -  02(1. K)  .  UCON  +  P  .  DT  •  0.25  •  (Z(I.K+1)  -  Z(I.K-1)) 

RHS(I.3)  -  03(1. K)  •  UCON  -  P  .  DT  .  0.25  •  (X( I .K+1 )-X( I ,K-1 ) ) 

RHS0.4)  =  UCON  •  (04(1.  K)+P)  -  XIT  »  P 

10  CONTINUE 

DO  1 1  1-2  .  IMAX  -  1 

D01(I.K)  =  D01(I.K)  -  RHS( 1+1,1)  +  RHS( 1-1.1) 
DQ2(I.K)  ■=  DQ2(I.K)  -  RHS(I  +  1.2)  +  RHS(I-1.2) 
DQ3(I.K)  =  DQ3(I.K)  -  RHS(I+1.3)  +  RHS(I-1.3) 

11  D04(I.K)  -  DQ4(I.K)  -  RHS(I+1.4)  +  RHS(I-1.4) 
100  CONTINUE 

C 

C       FLUX  TERMS  WITHIN  THE  ETA-  DERIVATIVE 

C 

DO  200  1=2.  IMAX  -  1 

DO  20  K  =  1  ,  KMAX 

VCON  -  (02(1 ,K)/01(I ,K)  )  •  (Z(I-1 ,K)-Z(I+1 .K)) 
1       +(Q3(I ,K)/Q1(I ,K)  )  •  (X(I+1 ,K)-X(I-1 ,K)) 

VCON  =  VCON  •  0.25  •  DT 

ETAT  =  -XTAU(I.K)  .  (Z( 1-1 ,K)-Z( 1+1 ,K) )  -  YTAU(I.K). 
1  (X(I+1 ,K)-X(I-1 ,K)) 

ETAT  =  ETAT  •  0.25  •  DT 

VCON  -  VCON  +  ETAT 

RHS(K. 1)  =  VCON  .01(1 ,K) 

P  =  (GAMMA-1.)  «  (04(1. K)  -  0.5  • (Q2( I ,K) ..2+03( 1 . K) ..2)/Q1 ( I . K)) 

RHS(K.2)  =  VCON  •  Q2(I,K)  +P  •  DT  .  .25  •  (Z( 1-1 ,K)-Z( 1+1 ,K) ) 

RHS(K,3)  =  VCON  •  03(1, K)  +  P  •  DT  .  . 25  •  (X(I+1,K)  -  X(I-I.K)) 

RHS(K,4)  =  VCON  •  (Q4( I. K)+P)  -  ETAT  «  P 

20  CONTINUE 

DO  21  K  =  2  .  KMAX  -  1 

D01(I.K)  -  D01(I.K)  -  RHS(K+1.1)  +  RHS(K-1.1) 
DQ2(I.K)  -  DQ2(I.K)  -  RHS(K+1 ,2)  +  RHS(K-1 .2) 
DQ3(I.K)  =  D03(I.K)  -  RHS(K+1 .3)  +  RHS(K-1,3) 

21  D04(I.K)  -  DQ4(I.K)  -  RHS(K+1.4)  +  RHS(K-1,4) 
200  CONTINUE 

RETURN 

END 
C 
C. ........................*.................... ......................... 

c 

SUBROUTINE  ROTGRID(X.Z. IMAX .KMAX .DALFA) 
C      ROTATE  GRID  IN  THE  CLOCKWISE  DIRECTION  BY  AN  AMOUNT  DALFA 
DIMENSION  X(161 .41 ) .Z( 161 .41) 
CA  =  COS (DALFA) 
SA  =  -  SIN(DALFA) 
DO  20  K  «  1  .  KMAX 
DO  20  I  -  1  .  IMAX 
X1  «=  X(I  ,K) 
Z1  =  Z(I.K) 
X(I  ,K)  «  X1  •  CA  -  Z1  •  SA 


20  Z(I ,K)  «  Z1  *  CA  +  X1  »  SA 
RETURN 
END 
C 
C. 

c 

SUBROUTINE  CPPL0T(I1 , 12 . FMACH.X, Y.CP) 
C 

C      THIS  SUBROUTINE  PLOTS  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE 
C 
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COMMON/SK  INCr/CF(  161) 

DIMENSION  KODE(4).LINE(90).X(161).Y(161) ,CP(161) 
DIMENSION  CFX(3.49).CFY(3.49),CPX(3.49).CPY(3.49) 
DATA  KODE/1H  .1H+.1HI ,1H«/ 
.      WRITE  (6.2) 

•  2  FORMA T(50H0P LOT  OF  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE/ 

•  1        10H0     X/C   .10H      CPU   ,10H      CFU   . 1 0H      CFL   ) 
CP0  -  (1.  +  .2  •  FMACH  ««2)  ••  3.5  -  1. 

CP0  -  CP0  /  (  .7  •  FMACH  «»2) 
K0  -  30.  •  CP0  +4.5 
IMIN  -  (I2-H)/2  +  11 
ILOW  =  2  •  IMIN 
I  COUNT  -  0 
CHD-X(H)  -  X(IMIN) 
DO  12  I  -  1  .90 
12  LINE(I)  -  KODE(1) 
DO  34  I  -  IMIN   .  12 
K  =  30.  •  (CP0  -  CP(I))  +  4.5 
K1  =  30.  •  (CP0  -CP(ILOW-I))  +  4.5 
IF(K.LT.I)  K  -  1 
IF(K1 . LT . 1 )  K1  -  1 
IF(K.GT.90)  K  -  90 
IF(K1  .GT.  90)  K1  -  90 
LINE(K0)  -  KODE(3) 
LINE(K)  -  KODE(2) 
LINE(K1)  -  K0DE(4) 
XOC  -  (X(I)  -  X(IMIN))  /  CHD 
•      WRITE  (6.610)  XOC. CP(I).CF(I),CF(ILOW-I). LINE 
LINE(K1)  -  KODE(1) 
34  LlNE(K)  -  KODE(1) 
C»»»   GENERATE  PLOT3D  CP   AND  CF  PLOTTING  FILES 
DO  500  I-IMIN.I2 

XOC  =  (X(I)  -  X(IMIN))/CHD 
I  COUNT  -  I COUNT  +  1 
CPY(1 . I COUNT)  -  0.000000 
CFY(1 , I COUNT)  -  0.000000 
CPY(2.ICOUNT)  -  CP(ILOW  -  I) 
CFY(2.ICOUNT)  -  CF(ILOW-1) 
CPY(3.ICOUNT)  -  CP(1) 
CFY(3.ICOUNT)  -  CF(I) 
CPX(1 .ICOUNT)  -  XOC 
CFX(1 . ICOUNT)  «  XOC 
CPX(2. ICOUNT)  =  XOC 
CFX(2. ICOUNT)  -  XOC 
CPX(3. ICOUNT)  -  XOC 
CFX(3. ICOUNT)  -  XOC 
500    CONTINUE 
IDM  -=  3 

JDM  =  12  -  IMIN  +  1 
WRITE(50)  IDM. JDM 

WRITE(50)((CPX(I,J).  1-1 .IDM) .J-1 .JDM). 
+  ((CPY(I.J),  1-1. IDM). J-1 .JDM) 

WRITEr60)  IDM. JDM 

WRITE(60)((CFX( I.J).  1-1 . IDM), J-1 .JDM). 
+  ((CFY(I.J).  1-1, IDM). J-1 .JDM) 

RETURN 
610  FORMAT(4F10.4,90A1) 
END 
/EOF 
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TITLE: 

NACA  0012  AIRFOIL 
I  MAX:   KMAX:   DT :      WW: 
161     41      .005     5. 
ITEL:   ITEU:   REYREF :  DNMIN: 
31      127    3  45     00005 
CSTP:   CPLT:   NSTP :   PSTP: 
0      0       1000    50 
FNU:       FNL:       FSYM: 
33.       33.       1. 
X:         Y: 


ALFA:  ALFA1 :  ALFAI:   REDFRE:  AMINF 

15.00  10.00  19.00   0.151    .283 

TSTART:  FORMAT:  RSTRT:  PITCH:  PLUNGE: 

-1.0  3.0  TRUE  TRUE    FALSE 


i.  0050 
.0125 
.0250 
.0500 
.0750 
.  1000 
.  1560 
.2000 
.2500 
.2600 
.2700 
.2800 
.2900 
.3000 
.3100 
.3200 
.3300 
.3400 
.3500 
.4000 
.4500 
.5000 
.5500 
.6000 
.6500 
.7000 
.7500 
.8000 
.8500 
.9000 
.9500 
.0000 


0 


.01153 
.01894 
.02615 
.03555 
.04200 
.04683 
.05345 
.05737 
.05941 
.05962 
.05978 
.05992 
.05999 
.06000 
.05999 
.05992 
.05980 
05965 
.05947 
.05803 
.05581 
.05294 
.04952 
.04563 
.04137 
.03664 
.03161 
.02623 
.02055 
.01448 
.00807 
.00126 
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